Data_analysis_MengHsinMorrisWu
Meng-Hsin, Wu
2023-07-24
1.1. Data manipulation on regional composition (Kenting)
For each site, three types of data will be arranged. The first type is a composition data based on initial categories; the second type is a composition data based on simplified categories; the third type is a data of substrate coverages.
#Get compositions in Kenting (KT)
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/Kenting")
library(readxl)
#Import data of Tiaoshih
TSdata<-read_xlsx('TS_composition.xlsx') #Biotic benthic cover
TSdata<-TSdata[,c(3,4,5,7,8)]
TSdata2<-read_xlsx('TS_composition.xlsx', sheet=2) #Substrate cover
TSdata2<-TSdata2[,c(3,4,6,8)]
TSdata<-as.data.frame(TSdata)
TSdata2<-as.data.frame(TSdata2)
TSdata.sum<-aggregate(TSdata[,3],by=list(Mfgroup=TSdata$Mfgroup), FUN=sum )
TSdata.sum.MFG2<-aggregate(TSdata[,3],
by=list(Mfgroup=TSdata$MFG2), FUN=sum )
prop<-TSdata.sum[,2]/TSdata2[1,3] #Get coverage
prop2<-TSdata.sum.MFG2[,2]/TSdata2[1,3]#Get coverage
TSdata.sum<-cbind(TSdata.sum,prop)
TSdata.sum.MFG2<-cbind(TSdata.sum.MFG2,prop2)
TSdata.sum[4]<-rep('TS', nrow(TSdata.sum))
TSdata.sum.MFG2[4]<-rep('TS', nrow(TSdata.sum.MFG2))
colnames(TSdata.sum)[4]<-'site'
colnames(TSdata.sum.MFG2)[4]<-'site'
write.csv(TSdata.sum, 'TS_sum.csv') #save table
write.csv(TSdata.sum.MFG2, 'TS_sum_MFG2.csv') #save table
TS2.sum<-aggregate(TSdata2[,3], by=list(Mfgroup=TSdata2$Mfgroup), FUN=sum )#Aggregate substrate coverage
prop<-TS2.sum[,2]/TSdata2[1,3]
TS2.sum$prop<-prop
TS2.sum$site<-rep('TS', nrow(TS2.sum))
write.csv(TS2.sum,'TS_substrates_sum.csv')#Save substrate coverage
###################################################
#Import data of Jialeshui
JLSdata<-read_xlsx('JLS_composition.xlsx')
JLSdata<-JLSdata[,c(2,3,4,6,7)]
JLSdata2<-read_xlsx('JLS_composition.xlsx', sheet=2)
JLSdata2<-JLSdata2[,c(3,4,6,8)]
JLSdata<-as.data.frame(JLSdata)
JLSdata2<-as.data.frame(JLSdata2)
JLSdata.sum<-aggregate(JLSdata[,3],by=list(Mfgroup=JLSdata$Mfgroup), FUN=sum )
JLSdata.sum.MFG2<-aggregate(JLSdata[,3],by=list(Mfgroup=JLSdata$MFG2), FUN=sum )
prop<-JLSdata.sum[,2]/JLSdata2[1,3]
prop2<-JLSdata.sum.MFG2[,2]/JLSdata2[1,3]
JLSdata.sum<-cbind(JLSdata.sum, prop)
JLSdata.sum.MFG2<-cbind(JLSdata.sum.MFG2, prop2)
JLSdata.sum[4]<-rep('JLS', nrow(JLSdata.sum))
JLSdata.sum.MFG2[4]<-rep('JLS', nrow(JLSdata.sum.MFG2))
colnames(JLSdata.sum)[4]<-'site'
colnames(JLSdata.sum.MFG2)[4]<-'site'
write.csv(JLSdata.sum, 'JLS_sum.csv') #save table
write.csv(JLSdata.sum.MFG2, 'JLS_sum_MFG2.csv') #save table
#Aggregate substrate coverage
JLS2.sum<-aggregate(JLSdata2[,3],
by=list(Mfgroup=JLSdata2$Mfgroup), FUN=sum )
prop<-JLS2.sum[,2]/JLSdata2[1,3]
JLS2.sum$prop<-prop
JLS2.sum$site<-rep('JLS', nrow(JLS2.sum))
write.csv(JLS2.sum,'JLS_substrates_sum.csv')
###############################################
#Import data of Houbihu
HBHdata<-read_xlsx('HBH_composition.xlsx')
HBHdata<-HBHdata[,c(2,3,4,6,7)]
HBHdata2<-read_xlsx('HBH_composition.xlsx', sheet=2)
HBHdata2<-HBHdata2[,c(3,4,6,8)]
HBHdata<-as.data.frame(HBHdata)
HBHdata2<-as.data.frame(HBHdata2)
HBHdata.sum<-aggregate(HBHdata[,3],by=list(Mfgroup=HBHdata$Mfgroup), FUN=sum )
HBHdata.sum.MFG2<-aggregate(HBHdata[,3],by=list(Mfgroup=HBHdata$MFG2), FUN=sum )
prop<-HBHdata.sum[,2]/HBHdata2[1,3]
prop2<-HBHdata.sum.MFG2[,2]/HBHdata2[1,3]
HBHdata.sum<-cbind(HBHdata.sum, prop)
HBHdata.sum.MFG2<-cbind(HBHdata.sum.MFG2, prop2)
HBHdata.sum<-as.data.frame(HBHdata.sum)
HBHdata.sum[4]<-rep('HBH', nrow(HBHdata.sum))
HBHdata.sum.MFG2[4]<-rep('HBH', nrow(HBHdata.sum.MFG2))
colnames(HBHdata.sum)[4]<-'site'
colnames(HBHdata.sum.MFG2)[4]<-'site'
write.csv(HBHdata.sum, 'HBH_sum.csv') #save table
write.csv(HBHdata.sum.MFG2, 'HBH_sum_MFG2.csv') #save table
#Aggregate substrate coverage
HBH2.sum<-aggregate(HBHdata2[,3],
by=list(Mfgroup=HBHdata2$Mfgroup), FUN=sum )
prop<-HBH2.sum[,2]/HBHdata2[1,3]
HBH2.sum$prop<-prop
HBH2.sum$site<-rep('HBH', nrow(HBH2.sum))
write.csv(HBH2.sum,'HBH_substrates_sum.csv')
################################################
#Import data of Leidashih
LDSdata<-read_xlsx('LDS_composition.xlsx')
LDSdata<-LDSdata[,c(3,4,5,7,8)]
LDSdata2<-read_xlsx('LDS_composition.xlsx', sheet=2)
LDSdata2<-LDSdata2[,c(3,4,6)]
LDSdata<-as.data.frame(LDSdata)
LDSdata2<-as.data.frame(LDSdata2)
LDSdata.sum<-aggregate(LDSdata[,3],
by=list(Mfgroup=LDSdata$Mfgroup), FUN=sum )
LDSdata.sum.MFG2<-aggregate(LDSdata[,3],
by=list(Mfgroup=LDSdata$MFG2), FUN=sum )
prop<-LDSdata.sum[,2]/LDSdata2[1,2]
prop2<-LDSdata.sum.MFG2[,2]/LDSdata2[1,2]
LDSdata.sum<-cbind(LDSdata.sum,prop)
LDSdata.sum.MFG2<-cbind(LDSdata.sum.MFG2,prop2)
LDSdata.sum[4]<-rep('LDS', nrow(LDSdata.sum))
LDSdata.sum.MFG2[4]<-rep('LDS', nrow(LDSdata.sum.MFG2))
colnames(LDSdata.sum)[4]<-'site'
colnames(LDSdata.sum.MFG2)[4]<-'site'
write.csv(LDSdata.sum, 'LDS_sum.csv') #save table
write.csv(LDSdata.sum.MFG2, 'LDS_sum_MFG2.csv') #save table
#Aggregate substrate coverage
LDS2.sum<-aggregate(LDSdata2[,2],
by=list(Mfgroup=LDSdata2$Mfgroup), FUN=sum )
prop<-LDS2.sum[,2]/LDSdata2[1,2]
LDS2.sum$prop<-prop
LDS2.sum$site<-rep('LDS', nrow(LDS2.sum))
write.csv(LDS2.sum,'LDS_substrates_sum.csv')
####################################################
#Import data of Dinbaisha
DiBSdata<-read_xlsx('DiBS_composition.xlsx')
DiBSdata<-DiBSdata[,c(3,4,5,7,8)]
DiBSdata2<-read_xlsx('DiBS_composition.xlsx', sheet=2)
DiBSdata2<-DiBSdata2[,c(3,4,6,8)]
DiBSdata<-as.data.frame(DiBSdata)
DiBSdata2<-as.data.frame(DiBSdata2)
DiBSdata.sum<-aggregate(DiBSdata[,3],
by=list(Mfgroup=DiBSdata$Mfgroup), FUN=sum )
DiBSdata.sum.MFG2<-aggregate(DiBSdata[,3],
by=list(Mfgroup=DiBSdata$MFG2), FUN=sum )
prop<-DiBSdata.sum[,2]/DiBSdata2[1,3]
prop2<-DiBSdata.sum.MFG2[,2]/DiBSdata2[1,3]
DiBSdata.sum<-cbind(DiBSdata.sum, prop)
DiBSdata.sum.MFG2<-cbind(DiBSdata.sum.MFG2, prop2)
DiBSdata.sum[4]<-rep('DiBS', nrow(DiBSdata.sum))
DiBSdata.sum.MFG2[4]<-rep('DiBS', nrow(DiBSdata.sum.MFG2))
colnames(DiBSdata.sum)[4]<-'site'
colnames(DiBSdata.sum.MFG2)[4]<-'site'
write.csv(DiBSdata.sum, 'DiBS_sum.csv') #save table
write.csv(DiBSdata.sum.MFG2, 'DiBS_sum_MFG2.csv') #save table
#Aggregate substrate coverage
DiBS2.sum<-aggregate(DiBSdata2[,3],
by=list(Mfgroup=DiBSdata2$Mfgroup), FUN=sum )
prop<-DiBS2.sum[,2]/DiBSdata2[1,3]
DiBS2.sum$prop<-prop
DiBS2.sum$site<-rep('DiBS', nrow(DiBS2.sum))
write.csv(DiBS2.sum,'DiBS_substrates_sum.csv')1.2. Data manipulation on regional composition (Lanyu)
#####Organizing LY data##############
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/Lanyu")
library(readxl)
###Import data of GreenGrassland
CCG<-read_xlsx('CCG_composition.xlsx')
CCG<-CCG[,c(3,4,6,8,9)]
CCG2<-read_xlsx('CCG_composition.xlsx', sheet=2)
CCG2<-CCG2[,c(3,4,6,8)]
CCG<-as.data.frame(CCG)
CCG2<-as.data.frame(CCG2)
CCG.sum<-aggregate(CCG[,3],
by=list(Mfgroup=CCG$Mfgroup), FUN=sum )
CCG.MFG2.sum<-aggregate(CCG[,3],
by=list(Mfgroup=CCG$MFG2), FUN=sum )
prop<-CCG.sum[,2]/CCG2[1,3]
prop<-as.data.frame(prop)
CCG.sum<-cbind(CCG.sum, prop)
CCG.sum[4]<-rep('CCG', nrow(CCG.sum))
colnames(CCG.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(CCG.sum,'CCG_sum.csv')
prop<-CCG.MFG2.sum[,2]/CCG2[1,3]
prop<-as.data.frame(prop)
CCG.MFG2.sum<-cbind(CCG.MFG2.sum, prop)
CCG.MFG2.sum[4]<-rep('CCG', nrow(CCG.MFG2.sum))
colnames(CCG.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(CCG.MFG2.sum,'CCG_MFG2_sum.csv')
#Substrate coverages
CCG2.sum<-aggregate(CCG2[,3],
by=list(Mfgroup=CCG2$Mfgroup), FUN=sum )
prop<-CCG2.sum[,2]/CCG2[1,3]
CCG2.sum$prop<-prop
CCG2.sum$site<-rep('CCG', nrow(CCG2.sum))
write.csv(CCG2.sum,'CCG_substrates_sum.csv')
######Import data of HenRock#
HR<-read_xlsx('HR_composition.xlsx')
HR<-HR[,c(3,4,6,8,9)]
HR2<-read_xlsx('HR_composition.xlsx', sheet=2)
HR2<-HR2[,c(3,4,6,8)]
HR<-as.data.frame(HR)
HR2<-as.data.frame(HR2)
HR.sum<-aggregate(HR[,3],
by=list(Mfgroup=HR$Mfgroup), FUN=sum )
HR.MFG2.sum<-aggregate(HR[,3],
by=list(Mfgroup=HR$MFG2), FUN=sum )
prop<-HR.sum[,2]/HR2[1,3]
prop<-as.data.frame(prop)
HR.sum<-cbind(HR.sum, prop)
HR.sum[4]<-rep('HR', nrow(HR.sum))
colnames(HR.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(HR.sum,'HR_sum.csv')
prop<-HR.MFG2.sum[,2]/HR2[1,3]
prop<-as.data.frame(prop)
HR.MFG2.sum<-cbind(HR.MFG2.sum, prop)
HR.MFG2.sum[4]<-rep('HR', nrow(HR.MFG2.sum))
colnames(HR.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(HR.MFG2.sum,'HR_MFG2_sum.csv')
#substrates coverage
HR2.sum<-aggregate(HR2[,3],
by=list(Mfgroup=HR2$Mfgroup), FUN=sum )
prop<-HR2.sum[,2]/HR2[1,3]
HR2.sum$prop<-prop
HR2.sum$site<-rep('HR', nrow(HR2.sum))
write.csv(HR2.sum,'HR_substrates_sum.csv')
##Import data of Dongchin#
DC<-read_xlsx('DC_composition.xlsx')
DC<-DC[,c(3,4,6,8,9)]
DC2<-read_xlsx('DC_composition.xlsx', sheet=2)
DC2<-DC2[,c(3,4,6,8)]
DC<-as.data.frame(DC)
DC2<-as.data.frame(DC2)
DC.sum<-aggregate(DC[,3],
by=list(Mfgroup=DC$Mfgroup), FUN=sum )
DC.MFG2.sum<-aggregate(DC[,3],
by=list(Mfgroup=DC$MFG2), FUN=sum )
prop<-DC.sum[,2]/DC2[1,3]
prop<-as.data.frame(prop)
DC.sum<-cbind(DC.sum, prop)
DC.sum[4]<-rep('DC', nrow(DC.sum))
colnames(DC.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(DC.sum,'DC_sum.csv')
prop<-DC.MFG2.sum[,2]/DC2[1,3]
prop<-as.data.frame(prop)
DC.MFG2.sum<-cbind(DC.MFG2.sum, prop)
DC.MFG2.sum[4]<-rep('DC', nrow(DC.MFG2.sum))
colnames(DC.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(DC.MFG2.sum,'DC_MFG2_sum.csv')
#substrate coverage#
DC2.sum<-aggregate(DC2[,3],
by=list(Mfgroup=DC2$Mfgroup), FUN=sum )
prop<-DC2.sum[,2]/DC2[1,3]
DC2.sum$prop<-prop
DC2.sum$site<-rep('DC', nrow(DC2.sum))
write.csv(DC2.sum,'DC_substrates_sum.csv')
##Import data of TwoLions#
TL<-read_xlsx('TL_composition.xlsx')
TL<-TL[,c(3,4,6,8,9)]
TL2<-read_xlsx('TL_composition.xlsx', sheet=2)
TL2<-TL2[,c(3,4,6,8)]
TL<-as.data.frame(TL)
TL2<-as.data.frame(TL2)
TL.sum<-aggregate(TL[,3],
by=list(Mfgroup=TL$Mfgroup), FUN=sum )
TL.MFG2.sum<-aggregate(TL[,3],
by=list(Mfgroup=TL$MFG2), FUN=sum )
prop<-TL.sum[,2]/TL2[1,3]
prop<-as.data.frame(prop)
TL.sum<-cbind(TL.sum, prop)
TL.sum[4]<-rep('TL', nrow(TL.sum))
colnames(TL.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(TL.sum,'TL_sum.csv')
prop<-TL.MFG2.sum[,2]/TL2[1,3]
prop<-as.data.frame(prop)
TL.MFG2.sum<-cbind(TL.MFG2.sum, prop)
TL.MFG2.sum[4]<-rep('TL', nrow(TL.MFG2.sum))
colnames(TL.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(TL.MFG2.sum,'TL_MFG2_sum.csv')
#substrates coverage#
TL2.sum<-aggregate(TL2[,3],
by=list(Mfgroup=TL2$Mfgroup), FUN=sum )
prop<-TL2.sum[,2]/TL2[1,3]
TL2.sum$prop<-prop
TL2.sum$site<-rep('TL', nrow(TL2.sum))
write.csv(TL2.sum,'TL_substrates_sum.csv')
##Import data of YaYo#
YY<-read_xlsx('YY_composition.xlsx')
YY<-YY[,c(3,4,6,8,9)]
YY2<-read_xlsx('YY_composition.xlsx', sheet=2)
YY2<-YY2[,c(3,4,6,8)]
YY<-as.data.frame(YY)
YY2<-as.data.frame(YY2)
YY.sum<-aggregate(YY[,3],
by=list(Mfgroup=YY$Mfgroup), FUN=sum )
YY.MFG2.sum<-aggregate(YY[,3],
by=list(Mfgroup=YY$MFG2), FUN=sum )
prop<-YY.sum[,2]/YY2[1,3]
prop<-as.data.frame(prop)
YY.sum<-cbind(YY.sum, prop)
YY.sum[4]<-rep('YY', nrow(YY.sum))
colnames(YY.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(YY.sum,'YY_sum.csv')
prop<-YY.MFG2.sum[,2]/YY2[1,3]
prop<-as.data.frame(prop)
YY.MFG2.sum<-cbind(YY.MFG2.sum, prop)
YY.MFG2.sum[4]<-rep('YY', nrow(YY.MFG2.sum))
colnames(YY.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(YY.MFG2.sum,'YY_MFG2_sum.csv')
#substrate coverage
YY2.sum<-aggregate(YY2[,3],
by=list(Mfgroup=YY2$Mfgroup), FUN=sum )
prop<-YY2.sum[,2]/YY2[1,3]
YY2.sum$prop<-prop
YY2.sum$site<-rep('YY', nrow(YY2.sum))
write.csv(YY2.sum,'YY_substrates_sum.csv')1.3. Data manipulation on regional composition (Ludao)
#####Organizing Ludao data##############
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/Ludao")
library(readxl)
##Import data of Guiwan
GW<-read_xlsx('GW_composition.xlsx')
GW<-GW[,c(3,4,6,8,9)]
GW2<-read_xlsx('GW_composition.xlsx', sheet=2)
GW2<-GW2[,c(3,4,6,8)]
GW<-as.data.frame(GW)
GW2<-as.data.frame(GW2)
GW.sum<-aggregate(GW[,3],
by=list(Mfgroup=GW$Mfgroup), FUN=sum )
GW.MFG2.sum<-aggregate(GW[,3],
by=list(Mfgroup=GW$MFG2), FUN=sum )
prop<-GW.sum[,2]/GW2[1,3]
prop<-as.data.frame(prop)
GW.sum<-cbind(GW.sum, prop)
GW.sum[4]<-rep('GW', nrow(GW.sum))
colnames(GW.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(GW.sum,'GW_sum.csv')
prop<-GW.MFG2.sum[,2]/GW2[1,3]
prop<-as.data.frame(prop)
GW.MFG2.sum<-cbind(GW.MFG2.sum, prop)
GW.MFG2.sum[4]<-rep('GW', nrow(GW.MFG2.sum))
colnames(GW.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(GW.MFG2.sum,'GW_MFG2_sum.csv')
#Substrate coverage
GW2.sum<-aggregate(GW2[,3],
by=list(Mfgroup=GW2$Mfgroup), FUN=sum )
prop<-GW2.sum[,2]/GW2[1,3]
GW2.sum$prop<-prop
GW2.sum$site<-rep('GW', nrow(GW2.sum))
write.csv(GW2.sum,'GW_substrates_sum.csv')
##Import data of Chaikou
CK<-read_xlsx('CK_composition.xlsx')
CK<-CK[,c(3,4,6,8,9)]
CK2<-read_xlsx('CK_composition.xlsx', sheet=2)
CK2<-CK2[,c(3,4,6,8)]
CK<-as.data.frame(CK)
CK2<-as.data.frame(CK2)
CK.sum<-aggregate(CK[,3],
by=list(Mfgroup=CK$Mfgroup), FUN=sum )
CK.MFG2.sum<-aggregate(CK[,3],
by=list(Mfgroup=CK$MFG2), FUN=sum )
prop<-CK.sum[,2]/CK2[1,3]
prop<-as.data.frame(prop)
CK.sum<-cbind(CK.sum, prop)
CK.sum[4]<-rep('CK', nrow(CK.sum))
colnames(CK.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(CK.sum,'CK_sum.csv')
prop<-CK.MFG2.sum[,2]/CK2[1,3]
prop<-as.data.frame(prop)
CK.MFG2.sum<-cbind(CK.MFG2.sum, prop)
CK.MFG2.sum[4]<-rep('CK', nrow(CK.MFG2.sum))
colnames(CK.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(CK.MFG2.sum,'CK_MFG2_sum.csv')
#substrate coverage
CK2.sum<-aggregate(CK2[,3],
by=list(Mfgroup=CK2$Mfgroup), FUN=sum )
prop<-CK2.sum[,2]/CK2[1,3]
CK2.sum$prop<-prop
CK2.sum$site<-rep('CK', nrow(CK2.sum))
write.csv(CK2.sum,'CK_substrates_sum.csv')
##Import data of Shihlang#
SL<-read_xlsx('SL_composition.xlsx')
SL<-SL[,c(3,4,6,8,9)]
SL2<-read_xlsx('SL_composition.xlsx', sheet=2)
SL2<-SL2[,c(3,4,6,8)]
SL<-as.data.frame(SL)
SL2<-as.data.frame(SL2)
SL.sum<-aggregate(SL[,3],
by=list(Mfgroup=SL$Mfgroup), FUN=sum )
SL.MFG2.sum<-aggregate(SL[,3],
by=list(Mfgroup=SL$MFG2), FUN=sum )
prop<-SL.sum[,2]/SL2[1,3]
prop<-as.data.frame(prop)
SL.sum<-cbind(SL.sum, prop)
SL.sum[4]<-rep('SL', nrow(SL.sum))
colnames(SL.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(SL.sum,'SL_sum.csv')
prop<-SL.MFG2.sum[,2]/SL2[1,3]
prop<-as.data.frame(prop)
SL.MFG2.sum<-cbind(SL.MFG2.sum, prop)
SL.MFG2.sum[4]<-rep('SL', nrow(SL.MFG2.sum))
colnames(SL.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(SL.MFG2.sum,'SL_MFG2_sum.csv')
#Substrate coverage
SL2.sum<-aggregate(SL2[,3],
by=list(Mfgroup=SL2$Mfgroup), FUN=sum )
prop<-SL2.sum[,2]/SL2[1,3]
SL2.sum$prop<-prop
SL2.sum$site<-rep('SL', nrow(SL2.sum))
write.csv(SL2.sum,'SL_substrates_sum.csv')
##Import Dabaisha data
DaBS<-read_xlsx('DaBS_composition.xlsx')
DaBS<-DaBS[,c(3,4,6,8,9)]
DaBS2<-read_xlsx('DaBS_composition.xlsx', sheet=2)
DaBS2<-DaBS2[,c(3,4,6,8)]
DaBS<-as.data.frame(DaBS)
DaBS2<-as.data.frame(DaBS2)
DaBS.sum<-aggregate(DaBS[,3],
by=list(Mfgroup=DaBS$Mfgroup), FUN=sum )
DaBS.MFG2.sum<-aggregate(DaBS[,3],
by=list(Mfgroup=DaBS$MFG2), FUN=sum )
prop<-DaBS.sum[,2]/DaBS2[1,3]
prop<-as.data.frame(prop)
DaBS.sum<-cbind(DaBS.sum, prop)
DaBS.sum[4]<-rep('DaBS', nrow(DaBS.sum))
colnames(DaBS.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(DaBS.sum,'DaBS_sum.csv')
prop<-DaBS.MFG2.sum[,2]/DaBS2[1,3]
prop<-as.data.frame(prop)
DaBS.MFG2.sum<-cbind(DaBS.MFG2.sum, prop)
DaBS.MFG2.sum[4]<-rep('DaBS', nrow(DaBS.MFG2.sum))
colnames(DaBS.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(DaBS.MFG2.sum,'DaBS_MFG2_sum.csv')
#Substrate coverage
DaBS2.sum<-aggregate(DaBS2[,3],
by=list(Mfgroup=DaBS2$Mfgroup), FUN=sum )
prop<-DaBS2.sum[,2]/DaBS2[1,3]
DaBS2.sum$prop<-prop
DaBS2.sum$site<-rep('DaBS', nrow(DaBS2.sum))
write.csv(DaBS2.sum,'DaBS_substrates_sum.csv')
##Import data of Gongguanbi#
GGB<-read_xlsx('GGB_composition.xlsx')
GGB<-GGB[,c(3,4,6,8,9)]
GGB2<-read_xlsx('GGB_composition.xlsx', sheet=2)
GGB2<-GGB2[,c(3,4,6,8)]
GGB<-as.data.frame(GGB)
GGB2<-as.data.frame(GGB2)
GGB.sum<-aggregate(GGB[,3],
by=list(Mfgroup=GGB$Mfgroup), FUN=sum )
GGB.MFG2.sum<-aggregate(GGB[,3],
by=list(Mfgroup=GGB$MFG2), FUN=sum )
prop<-GGB.sum[,2]/GGB2[1,3]
prop<-as.data.frame(prop)
GGB.sum<-cbind(GGB.sum, prop)
GGB.sum[4]<-rep('GGB', nrow(GGB.sum))
colnames(GGB.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(GGB.sum,'GGB_sum.csv')
prop<-GGB.MFG2.sum[,2]/GGB2[1,3]
prop<-as.data.frame(prop)
GGB.MFG2.sum<-cbind(GGB.MFG2.sum, prop)
GGB.MFG2.sum[4]<-rep('GGB', nrow(GGB.MFG2.sum))
colnames(GGB.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(GGB.MFG2.sum,'GGB_MFG2_sum.csv')
#Substrate coverage
GGB2.sum<-aggregate(GGB2[,3],
by=list(Mfgroup=GGB2$Mfgroup), FUN=sum )
prop<-GGB2.sum[,2]/GGB2[1,3]
GGB2.sum$prop<-prop
GGB2.sum$site<-rep('GGB', nrow(GGB2.sum))
write.csv(GGB2.sum,'GGB_substrates_sum.csv')1.4. Data manipulation on regional composition (East Coast)
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/East_Coast")
library(readxl)
##Import data of Jihui
JH<-read_xlsx('JH_composition.xlsx')
JH<-JH[,c(3,4,6,8,9)]
JH2<-read_xlsx('JH_composition.xlsx', sheet=2)
JH2<-JH2[,c(3,4,6,8)]
JH<-as.data.frame(JH)
JH2<-as.data.frame(JH2)
JH.sum<-aggregate(JH[,3],
by=list(Mfgroup=JH$Mfgroup), FUN=sum )
JH.MFG2.sum<-aggregate(JH[,3],
by=list(Mfgroup=JH$MFG2), FUN=sum )
prop<-JH.sum[,2]/JH2[1,3]
prop<-as.data.frame(prop)
JH.sum<-cbind(JH.sum, prop)
JH.sum[4]<-rep('JH', nrow(JH.sum))
colnames(JH.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(JH.sum,'JH_sum.csv')
prop<-JH.MFG2.sum[,2]/JH2[1,3]
prop<-as.data.frame(prop)
JH.MFG2.sum<-cbind(JH.MFG2.sum, prop)
JH.MFG2.sum[4]<-rep('JH', nrow(JH.MFG2.sum))
colnames(JH.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(JH.MFG2.sum,'JH_MFG2_sum.csv')
#Substrate coverage
JH2.sum<-aggregate(JH2[,3],
by=list(Mfgroup=JH2$Mfgroup), FUN=sum )
prop<-JH2.sum[,2]/JH2[1,3]
JH2.sum$prop<-prop
JH2.sum$site<-rep('JH', nrow(JH2.sum))
write.csv(JH2.sum,'JH_substrates_sum.csv')
##Import data of Jiamuziwan#
JMZ<-read_xlsx('JMZ_composition.xlsx')
JMZ<-JMZ[,c(3,4,6,8,9)]
JMZ2<-read_xlsx('JMZ_composition.xlsx', sheet=2)
JMZ2<-JMZ2[,c(3,4,6,8)]
JMZ<-as.data.frame(JMZ)
JMZ2<-as.data.frame(JMZ2)
JMZ.sum<-aggregate(JMZ[,3],
by=list(Mfgroup=JMZ$Mfgroup), FUN=sum )
JMZ.MFG2.sum<-aggregate(JMZ[,3],
by=list(Mfgroup=JMZ$MFG2), FUN=sum )
prop<-JMZ.sum[,2]/JMZ2[1,3]
prop<-as.data.frame(prop)
JMZ.sum<-cbind(JMZ.sum, prop)
JMZ.sum[4]<-rep('JMZ', nrow(JMZ.sum))
colnames(JMZ.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(JMZ.sum,'JMZ_sum.csv')
prop<-JMZ.MFG2.sum[,2]/JMZ2[1,3]
prop<-as.data.frame(prop)
JMZ.MFG2.sum<-cbind(JMZ.MFG2.sum, prop)
JMZ.MFG2.sum[4]<-rep('JMZ', nrow(JMZ.MFG2.sum))
colnames(JMZ.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(JMZ.MFG2.sum,'JMZ_MFG2_sum.csv')
#substrate coverage
JMZ2.sum<-aggregate(JMZ2[,3],
by=list(Mfgroup=JMZ2$Mfgroup), FUN=sum )
prop<-JMZ2.sum[,2]/JMZ2[1,3]
JMZ2.sum$prop<-prop
JMZ2.sum$site<-rep('JMZ', nrow(JMZ2.sum))
write.csv(JMZ2.sum,'JMZ_substrates_sum.csv')
##Import data of Fenniaolin#
FNL<-read_xlsx('FNL_composition.xlsx')
FNL<-FNL[,c(3,4,6,8,9)]
FNL2<-read_xlsx('FNL_composition.xlsx', sheet=2)
FNL2<-FNL2[,c(3,4,6,8)]
FNL<-as.data.frame(FNL)
FNL2<-as.data.frame(FNL2)
FNL.sum<-aggregate(FNL[,3],
by=list(Mfgroup=FNL$Mfgroup), FUN=sum )
FNL.MFG2.sum<-aggregate(FNL[,3],
by=list(Mfgroup=FNL$MFG2), FUN=sum )
prop<-FNL.sum[,2]/FNL2[1,3]
prop<-as.data.frame(prop)
FNL.sum<-cbind(FNL.sum, prop)
FNL.sum[4]<-rep('FNL', nrow(FNL.sum))
colnames(FNL.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(FNL.sum,'FNL_sum.csv')
prop<-FNL.MFG2.sum[,2]/FNL2[1,3]
prop<-as.data.frame(prop)
FNL.MFG2.sum<-cbind(FNL.MFG2.sum, prop)
FNL.MFG2.sum[4]<-rep('FNL', nrow(FNL.MFG2.sum))
colnames(FNL.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(FNL.MFG2.sum,'FNL_MFG2_sum.csv')
#substrate coverage
FNL2.sum<-aggregate(FNL2[,3],
by=list(Mfgroup=FNL2$Mfgroup), FUN=sum )
prop<-FNL2.sum[,2]/FNL2[1,3]
FNL2.sum$prop<-prop
FNL2.sum$site<-rep('FNL', nrow(FNL2.sum))
write.csv(FNL2.sum,'FNL_substrates_sum.csv')
##Import data of Hsinshe#
HS<-read_xlsx('HS_composition.xlsx')
HS<-HS[,c(3,4,6,8,9)]
HS2<-read_xlsx('HS_composition.xlsx', sheet=2)
HS2<-HS2[,c(3,4,6,8)]
HS<-as.data.frame(HS)
HS2<-as.data.frame(HS2)
HS.sum<-aggregate(HS[,3],
by=list(Mfgroup=HS$Mfgroup), FUN=sum )
HS.MFG2.sum<-aggregate(HS[,3],
by=list(Mfgroup=HS$MFG2), FUN=sum )
prop<-HS.sum[,2]/HS2[1,3]
prop<-as.data.frame(prop)
HS.sum<-cbind(HS.sum, prop)
HS.sum[4]<-rep('HS', nrow(HS.sum))
colnames(HS.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(HS.sum,'HS_sum.csv')
prop<-HS.MFG2.sum[,2]/HS2[1,3]
prop<-as.data.frame(prop)
HS.MFG2.sum<-cbind(HS.MFG2.sum, prop)
HS.MFG2.sum[4]<-rep('HS', nrow(HS.MFG2.sum))
colnames(HS.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(HS.MFG2.sum,'HS_MFG2_sum.csv')
#Substrate coverage
HS2.sum<-aggregate(HS2[,3],
by=list(Mfgroup=HS2$Mfgroup), FUN=sum )
prop<-HS2.sum[,2]/HS2[1,3]
HS2.sum$prop<-prop
HS2.sum$site<-rep('HS', nrow(HS2.sum))
write.csv(HS2.sum,'HS_substrates_sum.csv')
##Import data of Jiqi
JQ<-read_xlsx('JQ_composition.xlsx')
JQ<-JQ[,c(3,4,6,8,9)]
JQ2<-read_xlsx('JQ_composition.xlsx', sheet=2)
JQ2<-JQ2[,c(3,4,6,8)]
JQ<-as.data.frame(JQ)
JQ2<-as.data.frame(JQ2)
JQ.sum<-aggregate(JQ[,3],
by=list(Mfgroup=JQ$Mfgroup), FUN=sum )
JQ.MFG2.sum<-aggregate(JQ[,3],
by=list(Mfgroup=JQ$MFG2), FUN=sum )
prop<-JQ.sum[,2]/JQ2[1,3]
prop<-as.data.frame(prop)
JQ.sum<-cbind(JQ.sum, prop)
JQ.sum[4]<-rep('JQ', nrow(JQ.sum))
colnames(JQ.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(JQ.sum,'JQ_sum.csv')
prop<-JQ.MFG2.sum[,2]/JQ2[1,3]
prop<-as.data.frame(prop)
JQ.MFG2.sum<-cbind(JQ.MFG2.sum, prop)
JQ.MFG2.sum[4]<-rep('JQ', nrow(JQ.MFG2.sum))
colnames(JQ.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(JQ.MFG2.sum,'JQ_MFG2_sum.csv')
#substrate coverage
JQ2.sum<-aggregate(JQ2[,3],
by=list(Mfgroup=JQ2$Mfgroup), FUN=sum )
prop<-JQ2.sum[,2]/JQ2[1,3]
JQ2.sum$prop<-prop
JQ2.sum$site<-rep('JQ', nrow(JQ2.sum))
write.csv(JQ2.sum,'JQ_substrates_sum.csv')1.5. Data manipulation on regional composition (North Coast)
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment/North_Coast")
library(readxl)
##Import data of Bitou#
BT<-read_xlsx('BT_composition.xlsx')
BT<-BT[,c(3,4,6,8,9)]
BT2<-read_xlsx('BT_composition.xlsx', sheet=2)
BT2<-BT2[,c(3,4,6,8)]
BT<-as.data.frame(BT)
BT2<-as.data.frame(BT2)
BT.sum<-aggregate(BT[,3],
by=list(Mfgroup=BT$Mfgroup), FUN=sum )
BT.MFG2.sum<-aggregate(BT[,3],
by=list(Mfgroup=BT$MFG2), FUN=sum )
prop<-BT.sum[,2]/BT2[1,3]
prop<-as.data.frame(prop)
BT.sum<-cbind(BT.sum, prop)
BT.sum[4]<-rep('BT', nrow(BT.sum))
colnames(BT.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(BT.sum,'BT_sum.csv')
prop<-BT.MFG2.sum[,2]/BT2[1,3]
prop<-as.data.frame(prop)
BT.MFG2.sum<-cbind(BT.MFG2.sum, prop)
BT.MFG2.sum[4]<-rep('BT', nrow(BT.MFG2.sum))
colnames(BT.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(BT.MFG2.sum,'BT_MFG2_sum.csv')
#substrate coverage
BT2.sum<-aggregate(BT2[,3],
by=list(Mfgroup=BT2$Mfgroup), FUN=sum )
prop<-BT2.sum[,2]/BT2[1,3]
BT2.sum$prop<-prop
BT2.sum$site<-rep('BT', nrow(BT2.sum))
write.csv(BT2.sum,'BT_substrates_sum.csv')
##Import data of Chaojing#
CJ<-read_xlsx('CJ_composition.xlsx')
CJ<-CJ[,c(3,4,6,8,9)]
CJ2<-read_xlsx('CJ_composition.xlsx', sheet=2)
CJ2<-CJ2[,c(3,4,6,8)]
CJ<-as.data.frame(CJ)
CJ2<-as.data.frame(CJ2)
CJ.sum<-aggregate(CJ[,3],
by=list(Mfgroup=CJ$Mfgroup), FUN=sum )
CJ.MFG2.sum<-aggregate(CJ[,3],
by=list(Mfgroup=CJ$MFG2), FUN=sum )
prop<-CJ.sum[,2]/CJ2[1,3]
prop<-as.data.frame(prop)
CJ.sum<-cbind(CJ.sum, prop)
CJ.sum[4]<-rep('CJ', nrow(CJ.sum))
colnames(CJ.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(CJ.sum,'CJ_sum.csv')
prop<-CJ.MFG2.sum[,2]/CJ2[1,3]
prop<-as.data.frame(prop)
CJ.MFG2.sum<-cbind(CJ.MFG2.sum, prop)
CJ.MFG2.sum[4]<-rep('CJ', nrow(CJ.MFG2.sum))
colnames(CJ.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(CJ.MFG2.sum,'CJ_MFG2_sum.csv')
#substrate coverage
CJ2.sum<-aggregate(CJ2[,3],
by=list(Mfgroup=CJ2$Mfgroup), FUN=sum )
prop<-CJ2.sum[,2]/CJ2[1,3]
CJ2.sum$prop<-prop
CJ2.sum$site<-rep('CJ', nrow(CJ2.sum))
write.csv(CJ2.sum,'CJ_substrates_sum.csv')
##Import data of Longdong
LD<-read_xlsx('LD_composition.xlsx')
LD<-LD[,c(3,4,6,8,9)]
LD2<-read_xlsx('LD_composition.xlsx', sheet=2)
LD2<-LD2[,c(3,4,6,8)]
LD<-as.data.frame(LD)
LD2<-as.data.frame(LD2)
LD.sum<-aggregate(LD[,3],
by=list(Mfgroup=LD$Mfgroup), FUN=sum )
LD.MFG2.sum<-aggregate(LD[,3],
by=list(Mfgroup=LD$MFG2), FUN=sum )
prop<-LD.sum[,2]/LD2[1,3]
prop<-as.data.frame(prop)
LD.sum<-cbind(LD.sum, prop)
LD.sum[4]<-rep('LD', nrow(LD.sum))
colnames(LD.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(LD.sum,'LD_sum.csv')
prop<-LD.MFG2.sum[,2]/LD2[1,3]
prop<-as.data.frame(prop)
LD.MFG2.sum<-cbind(LD.MFG2.sum, prop)
LD.MFG2.sum[4]<-rep('LD', nrow(LD.MFG2.sum))
colnames(LD.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(LD.MFG2.sum,'LD_MFG2_sum.csv')
#substrate coverage
LD2.sum<-aggregate(LD2[,3],
by=list(Mfgroup=LD2$Mfgroup), FUN=sum )
prop<-LD2.sum[,2]/LD2[1,3]
LD2.sum$prop<-prop
LD2.sum$site<-rep('LD', nrow(LD2.sum))
write.csv(LD2.sum,'LD_substrates_sum.csv')
##Import data of 82.5K#
ZM<-read_xlsx('825_composition.xlsx')
ZM<-ZM[,c(3,4,6,8,9)]
ZM2<-read_xlsx('825_composition.xlsx', sheet=2)
ZM2<-ZM2[,c(3,4,6,8)]
ZM<-as.data.frame(ZM)
ZM2<-as.data.frame(ZM2)
ZM.sum<-aggregate(ZM[,3],
by=list(Mfgroup=ZM$Mfgroup), FUN=sum )
ZM.MFG2.sum<-aggregate(ZM[,3],
by=list(Mfgroup=ZM$MFG2), FUN=sum )
prop<-ZM.sum[,2]/ZM2[1,3]
prop<-as.data.frame(prop)
ZM.sum<-cbind(ZM.sum, prop)
ZM.sum[4]<-rep('82.5', nrow(ZM.sum))
colnames(ZM.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(ZM.sum,'82.5_sum.csv')
prop<-ZM.MFG2.sum[,2]/ZM2[1,3]
prop<-as.data.frame(prop)
ZM.MFG2.sum<-cbind(ZM.MFG2.sum, prop)
ZM.MFG2.sum[4]<-rep('82.5', nrow(ZM.MFG2.sum))
colnames(ZM.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(ZM.MFG2.sum,'82.5_MFG2_sum.csv')
#substrate coverage
ZM2.sum<-aggregate(ZM2[,3],
by=list(Mfgroup=ZM2$Mfgroup), FUN=sum )
prop<-ZM2.sum[,2]/ZM2[1,3]
ZM2.sum$prop<-prop
ZM2.sum$site<-rep('82.5', nrow(ZM2.sum))
write.csv(ZM2.sum,'82.5_substrates_sum.csv')
##Import data of Meiyenshan#
MYS<-read_xlsx('MYS_composition.xlsx')
MYS<-MYS[,c(3,4,6,8,9)]
MYS2<-read_xlsx('MYS_composition.xlsx', sheet=2)
MYS2<-MYS2[,c(3,4,6,8)]
MYS<-as.data.frame(MYS)
MYS2<-as.data.frame(MYS2)
MYS.sum<-aggregate(MYS[,3],
by=list(Mfgroup=MYS$Mfgroup), FUN=sum )
MYS.MFG2.sum<-aggregate(MYS[,3],
by=list(Mfgroup=MYS$MFG2), FUN=sum )
prop<-MYS.sum[,2]/MYS2[1,3]
prop<-as.data.frame(prop)
MYS.sum<-cbind(MYS.sum, prop)
MYS.sum[4]<-rep('MYS', nrow(MYS.sum))
colnames(MYS.sum)[c(2,4)]<-c('area','site')
#export table based on initial categories
write.csv(MYS.sum,'MYS_sum.csv')
prop<-MYS.MFG2.sum[,2]/MYS2[1,3]
prop<-as.data.frame(prop)
MYS.MFG2.sum<-cbind(MYS.MFG2.sum, prop)
MYS.MFG2.sum[4]<-rep('MYS', nrow(MYS.MFG2.sum))
colnames(MYS.MFG2.sum)[c(2,4)]<-c('area','site')
#export table based on simplified categories
write.csv(MYS.MFG2.sum,'MYS_MFG2_sum.csv')
#substrate coverage
MYS2.sum<-aggregate(MYS2[,3],
by=list(Mfgroup=MYS2$Mfgroup), FUN=sum )
prop<-MYS2.sum[,2]/MYS2[1,3]
MYS2.sum$prop<-prop
MYS2.sum$site<-rep('MYS', nrow(BT2.sum))
write.csv(MYS2.sum,'MYS_substrates_sum.csv')2. Data manipulation to combine all regional compositions.
This step is to combine all composition data with simplified categories, forming a big data frame. Within this data frame, each row tells one particular site, and each column tells one particular category.
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
######Import all data based on simplified categories
TS<-read.csv('~/myExperiment/Kenting/TS_sum_MFG2.csv')
LDS<-read.csv('~/myExperiment/Kenting/LDS_sum_MFG2.csv')
JLS<-read.csv('~/myExperiment/Kenting/JLS_sum_MFG2.csv')
HBH<-read.csv('~/myExperiment/Kenting/HBH_sum_MFG2.csv')
DiBS<-read.csv('~/myExperiment/Kenting/DiBS_sum_MFG2.csv')
GW<-read.csv('~/myExperiment/Ludao/GW_MFG2_sum.csv')
CK<-read.csv('~/myExperiment/Ludao/CK_MFG2_sum.csv')
SL<-read.csv('~/myExperiment/Ludao/SL_MFG2_sum.csv')
DaBS<-read.csv('~/myExperiment/Ludao/DaBS_MFG2_sum.csv')
GGB<-read.csv('~/myExperiment/Ludao/GGB_MFG2_sum.csv')
CCG<-read.csv('~/myExperiment/Lanyu/CCG_MFG2_sum.csv')
YY<-read.csv('~/myExperiment/Lanyu/YY_MFG2_sum.csv')
HR<-read.csv('~/myExperiment/Lanyu/HR_MFG2_sum.csv')
TL<-read.csv('~/myExperiment/Lanyu/TL_MFG2_sum.csv')
DC<-read.csv('~/myExperiment/Lanyu/DC_MFG2_Sum.csv')
JH<-read.csv('~/myExperiment/East Coast/JH_MFG2_sum.csv')
JMZ<-read.csv('~/myExperiment/East Coast/JMZ_MFG2_sum.csv')
HS<-read.csv('~/myExperiment/East Coast/HS_MFG2_sum.csv')
JQ<-read.csv('~/myExperiment/East Coast/JQ_MFG2_sum.csv')
FNL<-read.csv('~/myExperiment/East Coast/FNL_MFG2_sum.csv')
LD<-read.csv('~/myExperiment/North Coast/LD_MFG2_sum.csv')
BT<-read.csv('~/myExperiment/North Coast/BT_MFG2_sum.csv')
CJ<-read.csv('~/myExperiment/North Coast/CJ_MFG2_sum.csv')
MYS<-read.csv('~/myExperiment/North Coast/MYS_MFG2_sum.csv')
ZM<-read.csv('~/myExperiment/North Coast/82.5_MFG2_sum.csv')
###Organizing the datasets###
#Set the factor telling regions
colnames(TS)[3]<-c('area')
colnames(DiBS)[3]<-c('area')
colnames(JLS)[3]<-c('area')
colnames(LDS)[3]<-c('area')
colnames(HBH)[3]<-c('area')
colnames(TS)[4]<-c('prop')
colnames(DiBS)[4]<-c('prop')
colnames(JLS)[4]<-c('prop')
colnames(LDS)[4]<-c('prop')
colnames(HBH)[4]<-c('prop')
TS$region<-'KT'
DiBS$region<-'KT'
JLS$region<-'KT'
LDS$region<-'KT'
HBH$region<-'KT'
GW$region<-'LD'
CK$region<-'LD'
SL$region<-'LD'
DaBS$region<-'LD'
GGB$region<-'LD'
CCG$region<-'LY'
YY$region<-'LY'
HR$region<-'LY'
TL$region<-'LY'
DC$region<-'LY'
JH$region<-'EC'
JMZ$region<-'EC'
HS$region<-'EC'
JQ$region<-'EC'
FNL$region<-'EC'
LD$region<-'NC'
BT$region<-'NC'
CJ$region<-'NC'
MYS$region<-'NC'
ZM$region<-'NC'
#Get longer table
composition_long<-rbind(TS, LDS, JLS, HBH, DiBS,
GW,CK, SL, DaBS, GGB,
CCG, YY, HR, TL, DC,
JH, JMZ, HS, JQ, FNL,
LD, BT, CJ, MYS, ZM)
composition_long<-composition_long[c(2,4,5,6)]
#Get wider table
library(tidyr)
composition<-pivot_wider(composition_long, names_from = Mfgroup,
values_from = prop)
composition[is.na(composition)]<-0
composition<-as.data.frame(composition)
#Further simplify the categories
composition$other<-composition$other+composition$others+composition$crioth #combine crinoids with other invertebrates
composition$others<-NULL
composition$crioth<-NULL
composition$cru<-composition$cur+composition$cru #cur is a wrong category, it needs to be combined with cru
composition$cur<-NULL
composition$tws<-composition$wantws+composition$tws #combine all types of TWS together
composition$wantws<-NULL
#write.csv(composition, 'analysis_tables/Composition_AllSites.csv')
composition$cru_or_spoenc<-composition$cru+composition$spoenc #Since the difficulty to distinguish encrusting sponge and CCA, these two categories are merged as well.
composition$cru<-NULL
composition$spoenc<-NULL
write.csv(composition,'analysis_tables/Composition__CombineSpoencAndCru_AllSites.csv' )#####Importing substrate data########
TS2<-read.csv('~/myExperiment/Kenting/TS_substrates_sum.csv')
LDS2<-read.csv('~/myExperiment/Kenting/LDS_substrates_sum.csv')
JLS2<-read.csv('~/myExperiment/Kenting/JLS_substrates_sum.csv')
HBH2<-read.csv('~/myExperiment/Kenting/HBH_substrates_sum.csv')
DiBS2<-read.csv('~/myExperiment/Kenting/DiBS_substrates_sum.csv')
GW2<-read.csv('~/myExperiment/Ludao/GW_substrates_sum.csv')
CK2<-read.csv('~/myExperiment/Ludao/CK_substrates_sum.csv')
SL2<-read.csv('~/myExperiment/Ludao/SL_substrates_sum.csv')
DaBS2<-read.csv('~/myExperiment/Ludao/DaBS_substrates_sum.csv')
GGB2<-read.csv('~/myExperiment/Ludao/GGB_substrates_sum.csv')
CCG2<-read.csv('~/myExperiment/Lanyu/CCG_substrates_sum.csv')
YY2<-read.csv('~/myExperiment/Lanyu/YY_substrates_sum.csv')
HR2<-read.csv('~/myExperiment/Lanyu/HR_substrates_sum.csv')
TL2<-read.csv('~/myExperiment/Lanyu/TL_substrates_sum.csv')
DC2<-read.csv('~/myExperiment/Lanyu/DC_substrates_sum.csv')
JH2<-read.csv('~/myExperiment/East Coast/JH_substrates_sum.csv')
JMZ2<-read.csv('~/myExperiment/East Coast/JMZ_substrates_sum.csv')
HS2<-read.csv('~/myExperiment/East Coast/HS_substrates_sum.csv')
JQ2<-read.csv('~/myExperiment/East Coast/JQ_substrates_sum.csv')
FNL2<-read.csv('~/myExperiment/East Coast/FNL_substrates_sum.csv')
LD2<-read.csv('~/myExperiment/North Coast/LD_substrates_sum.csv')
BT2<-read.csv('~/myExperiment/North Coast/BT_substrates_sum.csv')
CJ2<-read.csv('~/myExperiment/North Coast/CJ_substrates_sum.csv')
MYS2<-read.csv('~/myExperiment/North Coast/MYS_substrates_sum.csv')
ZM2<-read.csv('~/myExperiment/North Coast/82.5_substrates_sum.csv')
#Make longer table
substrates_long<-rbind(TS2, LDS2, JLS2, HBH2, DiBS2,
GW2,CK2, SL2, DaBS2, GGB2,
CCG2, YY2, HR2, TL2, DC2,
JH2, JMZ2, HS2, JQ2, FNL2,
LD2, BT2, CJ2, MYS2, ZM2)
substrates_long<-substrates_long[c(2,4,5)]
#Make wider table
substrates<-pivot_wider(substrates_long, names_from = Mfgroup,
values_from = prop)
substrates[is.na(substrates)]<-0
substrates<-as.data.frame(substrates)
substrates$region<-composition$region
substrates<-substrates[,-2]
#Save the preliminary composition table
write.csv(substrates,'analysis_tables/Substrates_AllSites.csv' )3.1. Metric extraction (Fractal dimension)
This code is used to compute fractal dimension of different scale intervals based on surface area and planar area extracted from QGIS.
rm(list=ls())
setwd("~your working directory")
SA.data<-data.frame()
resolution<-c('The whole resolution', 1,2,4,8,16,32,64,128)
#Copy and paste surface area of different resolutions (accessed by function r.surf.area in GDAL)
SurfaceArea<-c( 'SA (whole Res)', 'SA (1 cm)', 'SA (2 cm)',
'SA (4 cm)', 'SA (8 cm)',
'SA (16 cm)', 'SA (32 cm)',
'SA (64 cm)', 'SA (128 cm)')
#Copy and paste calculated planar area of different resolutions.
PlanAreaCalc<-c( 'PAc (Whole Res)',
'PAc (1 cm)', 'PAc (2 cm)',
'PAc (4 cm)', 'PAc (8 cm)',
'PAc (16 cm)', 'PAc (32 cm)',
'PAc (64 cm)', ' PAc (128 cm)')
#Copy and paste null area (the true planar area is computed by subtracting calculated planar area from null area)
NullArea<-c( 'NA (Whole Res)',
'NA (1 cm)', 'NA (2 cm)',
'NA (4 cm)', 'NA (8 cm)',
'NA (16 cm)', 'NA (32 cm)',
'NA (64 cm)', ' NA (128 cm)')
PlanArea<-PlanAreaCalc-NullArea #computing the true planar area
SA.data<-data.frame(resolution=resolution,
SurfaceArea=SurfaceArea,
PlanArea=PlanArea)
SA.data$SurfaceComplexity<-
SA.data$SurfaceArea/SA.data$PlanArea
SA.data$LogSA<-log(SurfaceArea)
SA.data$LogRES<-log(resolution/100)
lmSA64<-lm(LogSA[2:8]~LogRES[2:8], data=SA.data)
FD64<-2-lmSA64_2$coefficients[2] # Get D (1 - 64 cm)
library(ggplot2)
#Visualization of D (1 - 64 cm)
FD64Plot<-ggplot(SA.data[1:8, ], aes(x=LogRES,
y=LogSA))+
geom_point()+
geom_smooth(method='lm')+
xlab('log(Resolution)')+
ylab('log(Surface Area)')+
ggtitle(paste('name of sampling site', 'D64= ', FD64_2))
####################################
#Get FD at different scale intervals
slopes<-c()
FD_intervals<-c()
for(i in 2:7){
slopes[i-1]<-(SA.data$LogSA[i+1]-SA.data$LogSA[i])/(SA.data$LogRES[i+1]-SA.data$LogRES[i])
}
FD_intervals= 2-slopes
intervals<-c('D1_2','D2_4', 'D4_8', 'D8_16', 'D16_32', 'D32_64')
FDs<-data.frame(interv= intervals, FDs= FD_intervals)3.2. Import metric table and conducting stepwise VIF selection on complexity metrics.
The metrics were collected by using QGIS. Fractal dimension (D) was computed by data extracted by QGIS and the code provided at section 1. All metrics of each site were manually organized in a table.
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
library(corrplot)
library(car)
library(readxl)
library(olsrr)
#Read metric data
StrMetrics<-read.csv('Sites_MetricSummary.csv', header=T)
row.names(StrMetrics)<-StrMetrics[,1]
metrics<-StrMetrics[, 2:24]
#Import relief data
relief_data<-read_xlsx('VRandVRSD.xlsx', sheet=2)
relief_data<-as.data.frame(relief_data)
row.names(relief_data)<-relief_data[,1]
relief_data<-relief_data[-1]
metrics<-merge(metrics, relief_data[2], by='row.names')
row.names(metrics)<-metrics$Row.names
metrics<-metrics[-1] #Remove VR, keep VRSD (Sq)
write.csv(metrics, 'Site_MetricSummary_2.csv') #Save the table with all metrics
head(metrics)## S4 T4 V4 S32 T32 V32 PROC4
## 82.5 35.79695 0.10297550 0.03972595 34.08894 0.6829777 0.04316221 3.138004
## BT 29.06211 0.07379687 0.06354869 18.46841 0.3176737 0.01612316 4.993615
## CCG 31.98746 0.09315995 0.08002010 22.27737 0.3669287 0.01298073 4.518048
## CJ 37.43711 0.13200221 0.07057097 29.69820 0.6569097 0.08023318 3.669758
## CK 39.58948 0.13459857 0.05272271 28.06464 0.6684020 0.09568557 2.722089
## DaBS 22.30837 0.05852110 0.04697374 10.64565 0.2187278 0.01104392 3.920315
## PLC4 PROC32 PLC32 S16 T16 V16 PROC16
## 82.5 9.242265 0.3542356 1.360989 34.22099 0.3514723 0.03748765 0.4714557
## BT 18.462807 0.3869983 1.543360 20.79732 0.2175307 0.03928283 1.0199474
## CCG 14.839175 0.2802628 1.027368 23.48149 0.2233868 0.02969492 0.8610387
## CJ 12.282814 0.4881684 1.838733 33.84235 0.3906941 0.07143632 0.9621377
## CK 8.519395 0.6309607 1.913845 34.12422 0.4080009 0.09667555 1.1375918
## DaBS 15.539315 0.2659814 1.967384 13.80435 0.1368370 0.01919955 0.7120210
## PLC16 D64 SC D1_2 D2_4 D4_8 D8_16 D16_32
## 82.5 5.380102 2.082288 1.637182 2.059848 2.059313 2.069394 2.048714 2.075094
## BT 5.696730 2.142668 1.617211 2.105806 2.136347 2.147361 2.146462 2.134220
## CCG 3.184027 2.141438 1.848779 2.196623 2.187020 2.149519 2.157569 2.034047
## CJ 3.296366 2.134910 1.979722 2.105106 2.102707 2.125147 2.165645 2.097948
## CK 3.775512 2.170443 2.080996 2.125388 2.149982 2.158457 2.136556 2.170571
## DaBS 3.822742 2.120286 1.374338 2.102362 2.106039 2.112706 2.097096 2.125669
## D32_64 VRSD4
## 82.5 2.247949 0.6966732
## BT 2.187176 0.3594060
## CCG 2.140841 0.5719821
## CJ 2.238046 0.4688934
## CK 2.341126 0.3662118
## DaBS 2.214526 0.1714850metrics_cor<-cor(metrics) #Get correlation matrix composed of every metric pair.
####################Corrplot adjustment####################################
metrics_pub<-metrics
colnames(metrics_pub)<-c('S - 4 cm', 'TRI - 4 cm', 'VRM - 4cm', 'S - 32 cm',
'TRI - 32 cm', 'VRM - 32 cm', 'PROC - 4 cm', 'PLC - 4 cm',
'PROC - 32 cm', 'PLC - 32 cm', 'S - 16 cm', 'TRI - 16 cm',
'VRM - 16 cm', 'PROC - 16 cm', 'PLC - 16 cm', 'D (1 - 64 cm)',
'SC', 'D (1 - 2 cm)', 'D (2 - 4 cm)', 'D (4 - 8 cm)', 'D (8 - 16 cm)',
'D (16 - 32 cm)', 'D (32 - 64 cm)', 'Sq')
metrics_cor<-cor(metrics_pub)
#corrplot(metrics_cor, method='number',
# type='lower')
res<-cor.mtest(metrics_cor,
conf.level = .95)
#########################Displaying Corrplot#############################
corrplot(metrics_cor, method='number',
type='lower', p.mat = res$p,
sig.level = .05, tl.cex=0.5, cl.cex=0.5,
number.cex=0.5)
#######################################################
###############Use usdm for stepwise VIF selection on metrics##################################
library(usdm)
v1<-vifstep(metrics, th=10)
metrics_selected<-exclude(metrics, v1)
metrics_selected #See metrics retained from selection## V4 PLC4 PROC32 PLC32 PLC16 D1_2 D2_4
## 82.5 0.03972595 9.242265 0.3542356 1.3609894 5.380102 2.059848 2.059313
## BT 0.06354869 18.462807 0.3869983 1.5433599 5.696730 2.105806 2.136347
## CCG 0.08002010 14.839175 0.2802628 1.0273676 3.184027 2.196623 2.187020
## CJ 0.07057097 12.282814 0.4881684 1.8387325 3.296366 2.105106 2.102707
## CK 0.05272271 8.519395 0.6309607 1.9138448 3.775512 2.125388 2.149982
## DaBS 0.04697374 15.539315 0.2659814 1.9673836 3.822742 2.102362 2.106039
## DC 0.07477847 17.502046 0.4166059 2.3384137 5.065466 2.175885 2.184300
## DiBS 0.04879224 13.948038 0.3448847 1.5995921 3.319096 2.146635 2.144978
## FNL 0.09231118 14.027762 0.5720152 1.8826472 3.643584 2.173418 2.185573
## GGB 0.06445081 14.743830 0.4614089 2.2478442 4.286704 2.142029 2.141366
## GW 0.04763611 15.339935 0.2874124 3.1911692 4.278659 2.125547 2.134165
## HBH 0.06507814 14.385386 0.5483947 1.7884567 3.607074 2.179179 2.169388
## HR 0.09157567 14.738506 0.5034706 2.2121285 4.682016 2.186643 2.206191
## HS 0.06418384 11.639053 0.7444362 2.4561036 3.899267 2.116186 2.117397
## JH 0.06854386 12.926524 0.5524131 1.9177405 3.275951 2.123935 2.120539
## JLS 0.08996626 13.711100 0.5674756 2.0121217 3.471119 2.213680 2.187870
## JMZ 0.10180321 13.580160 0.6749426 2.6838572 4.000385 2.147930 2.145423
## JQ 0.05014038 12.100439 0.5099651 2.0976212 2.811581 2.123612 2.102982
## LD 0.09773878 15.983127 0.2414503 0.8735963 4.529776 2.114213 2.142355
## LDS 0.07129418 18.708382 0.4316027 2.6060655 3.972630 2.117886 2.169715
## MYS 0.03591826 21.839115 0.1743400 1.7492916 4.682573 2.108662 2.116179
## SL 0.08072169 10.510869 0.5701214 1.8626170 3.517852 2.255038 2.152829
## TL 0.04717792 14.914221 0.2674686 1.8916910 4.247381 2.136491 2.124769
## TS 0.08013721 14.234224 0.4679286 2.5597380 4.639967 2.134288 2.130212
## YY 0.05256534 17.972904 0.2521757 1.9701828 5.083402 2.136654 2.139941
## D8_16 D16_32 D32_64 VRSD4
## 82.5 2.048714 2.075094 2.247949 0.69667321
## BT 2.146462 2.134220 2.187176 0.35940599
## CCG 2.157569 2.034047 2.140841 0.57198211
## CJ 2.165645 2.097948 2.238046 0.46889336
## CK 2.136556 2.170571 2.341126 0.36621178
## DaBS 2.097096 2.125669 2.214526 0.17148500
## DC 2.144677 2.146597 2.195416 0.21552578
## DiBS 2.093201 2.055113 2.128465 0.27035985
## FNL 2.164326 2.171833 2.256766 0.32194592
## GGB 2.153132 2.114266 2.312947 0.35606082
## GW 2.092675 2.130449 2.179806 0.09572413
## HBH 2.094460 2.137286 2.292320 0.47868320
## HR 2.160204 2.106320 2.213140 0.29349338
## HS 2.185740 2.188789 2.338200 0.31157591
## JH 2.148578 2.156038 2.229615 0.49583007
## JLS 2.171962 2.132067 2.255832 0.51166963
## JMZ 2.197287 2.243912 2.308156 0.32193789
## JQ 2.138418 2.138883 2.238538 0.47133979
## LD 2.182984 2.173301 2.164168 0.49713277
## LDS 2.149087 2.092245 2.201565 0.17018665
## MYS 2.059969 2.003714 2.069643 0.10205495
## SL 2.106764 2.121427 2.303921 0.35146390
## TL 2.065154 2.127289 2.087826 0.26951771
## TS 2.128885 2.163362 2.132155 0.22094171
## YY 2.055983 2.084683 2.153331 0.173852144.1. Comparison of complexity metrics (Latitudinal correlation)
##Latitudinal correlation of metrics###
library(readxl)
#Import latitudes of sites
coordi<-read_xlsx('Photogrammetric_Plot_Coordinates.xlsx') ## New names:
## • `` -> `...6`coordi<-as.data.frame(coordi)
coordi<-coordi[,-c(4,5,6)]
row.names(coordi)<-coordi[,1]
#Merge latitudes and metrics
metrics_coordi<-merge(metrics, coordi,
by='row.names')
row.names(metrics_coordi)<-metrics_coordi[,1]
metrics_coordi<-metrics_coordi[-1]
metrics_coordi<-merge(metrics_coordi, StrMetrics[25],
by='row.names')
metrics_coordi$region<-factor(metrics_coordi$region,
levels=c('KT','LY','LD',
'EC','NC'),
ordered = T)
#V4
cor_v4_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$V4, method='pearson',
conf.level = 0.95)
#Plot the relationship with latitudes
library(ggplot2)
plot_v4_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=V4))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=22.5, y=0.037, size=5,label=paste('p = ',round(cor_v4_lat$p.value, 3)))+
annotate('text', x=22.5, y=0.043, size=5,label=paste('r = ',round(cor_v4_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('VRM - 4 cm')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_v4_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_v4_lat.jpeg', plot_v4_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#PROC32
cor_PROC32_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$PROC32, method='pearson',
conf.level = 0.95)
#Plot the relationship with latitudes
plot_PROC32_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=PROC32))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=24.5, y=0.135, size=5, label=paste('p = ',round(cor_PROC32_lat$p.value, 3)))+
annotate('text', x=24.5, y=0.175, size=5, label=paste('r = ',round(cor_PROC32_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('PROC - 32 cm')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_PROC32_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_PROC32_lat.jpeg', plot_PROC32_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#PLC4
cor_PLC4_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$PLC4, method='pearson',
conf.level = 0.95)
#Plot the relationship
plot_PLC4_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=PLC4))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=22.5, y=3, size=5,label=paste('p = ',round(cor_PLC4_lat$p.value, 3)))+
annotate('text', x=22.5, y=5, size=5, label=paste('r = ',round(cor_PLC4_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('PLC - 4 cm')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_PLC4_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_PLC4_lat.jpeg', plot_PLC4_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#PLC16
cor_PLC16_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$PLC16, method='pearson',
conf.level = 0.95)
#Plot relationship
plot_PLC16_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=PLC16))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=22.5, y=2.9,size=5, label=paste('p = ',round(cor_PLC16_lat$p.value, 3)))+
annotate('text', x=22.5, y=3.2, size=5,label=paste('r = ',round(cor_PLC16_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('PLC - 16 cm')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_PLC16_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_PLC16_lat.jpeg', plot_PLC16_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#PLC32
cor_PLC32_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$PLC32, method='pearson',
conf.level = 0.95)
#no significant correlation btw sc and latitude
plot_PLC32_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=PLC32))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=22.6, y=1, size=5,label=paste('p = ',round(cor_PLC32_lat$p.value, 3)))+
annotate('text', x=22.6, y=1.2, size=5,label=paste('r = ',round(cor_PLC32_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('PLC - 32 cm')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_PLC32_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('figures/plot_PLC32_lat.jpeg', plot_PLC32_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#D12
cor_D12_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D1_2, method='pearson',
conf.level = 0.95)
#Plot relationship
plot_D12_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D1_2))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=22.5, y=2.07,size=6.5, label=paste('p = ',round(cor_D12_lat$p.value, 3)))+
annotate('text', x=22.5, y=2.085, size=6.5,label=paste('r = ',round(cor_D12_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('D (1 - 2 cm)')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D12_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_D12_lat.jpeg', plot_D12_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#D24
cor_D24_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D2_4, method='pearson',
conf.level = 0.95)
#Plot relationship
plot_D24_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D2_4))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=22.5, y=2.07, size=6.5, label=paste('p = ',round(cor_D24_lat$p.value, 3)))+
annotate('text', x=22.5, y=2.085, size=6.5, label=paste('r = ',round(cor_D24_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('D (2 - 4 cm)')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D24_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_D24_lat.jpeg', plot_D24_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#D816
cor_D816_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D8_16, method='pearson',
conf.level = 0.95)
#Plot relationship
plot_D816_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D8_16))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=23, y=2.057, size=5,label=paste('p = ',round(cor_D816_lat$p.value, 3)))+
annotate('text', x=23, y=2.07, size=5,label=paste('r = ',round(cor_D816_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('D (8 - 16 cm)')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D816_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_D816_lat.jpeg', plot_D816_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#D1632
cor_D1632_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D16_32, method='pearson',
conf.level = 0.95)
#Plot relationship
plot_D1632_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D16_32))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=23, y=2.02,size=5, label=paste('p = ',round(cor_D1632_lat$p.value, 3)))+
annotate('text', x=23, y=2.04, size=5,label=paste('r = ',round(cor_D1632_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('D (16 - 32 cm)')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D1632_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_D1632_lat.jpeg', plot_D1632_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#D3264
cor_D3264_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$D32_64, method='pearson',
conf.level = 0.95)
#Plot relationship
plot_D3264_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=D32_64))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=23, y=2.07, size=5,label=paste('p = ',round(cor_D3264_lat$p.value, 3)))+
annotate('text', x=23, y=2.09, size=5,label=paste('r = ',round(cor_D3264_lat$estimate, 3)))+
xlab('Latitudes')+ ylab('D (32 - 64 cm)')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_D3264_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_D3264_lat.jpeg', plot_D3264_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'#Sq
cor_VRSD_lat<-cor.test(metrics_coordi$Lat, metrics_coordi$VRSD4, method='pearson',
conf.level = 0.95)
#Plot relationship
plot_VRSD_lat<-ggplot(data=metrics_coordi, aes(x=Lat, y=VRSD4))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
annotate('text', x=23.5, y=0.15, size=5,label=paste('p = ',round(cor_VRSD_lat$p.value, 3)))+
annotate('text', x=23.5, y=0.18, size=5,label=paste('r = ',round(cor_VRSD_lat$estimate, 3)))+
labs(y='Sq')+xlab('Latitudes')+
guides(color = guide_legend(title = "Region"))+
scale_color_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
plot_VRSD_lat## `geom_smooth()` using formula = 'y ~ x'
ggsave('Figures/plot_VRSD_lat.jpeg', plot_VRSD_lat,
width=15, height=15, unit='cm')## `geom_smooth()` using formula = 'y ~ x'4.2. Comparison of complexity metrics (Interregional differences)
##########################Interregional comparisons: Kruskal Wallis test##############################
V4_kruskal<-kruskal.test(V4~region, data=metrics_coordi)
V32_kruskal<-kruskal.test(V32~region, data=metrics_coordi)
PROC32_kruskal<-kruskal.test(PROC32~region, data=metrics_coordi)
PLC4_kruskal<-kruskal.test(PLC4~region, data=metrics_coordi)
PLC16_kruskal<-kruskal.test(PLC16~region, data=metrics_coordi)
PLC32_kruskal<-kruskal.test(PLC32~region, data=metrics_coordi)
D12_kruskal<-kruskal.test(D1_2~region, data=metrics_coordi)
D24_kruskal<-kruskal.test(D2_4~region, data=metrics_coordi)
D48_kruskal<-kruskal.test(D4_8~region, data=metrics_coordi)
D816_kruskal<-kruskal.test(D8_16~region, data=metrics_coordi)
D1632_kruskal<-kruskal.test(D16_32~region, data=metrics_coordi)
D3264_kruskal<-kruskal.test(D32_64~region, data=metrics_coordi)
VRSD_kruskal<-kruskal.test(VRSD4~region, data=metrics_coordi)
######################Pairwise Dunn test########################################
library(PMCMRplus)
PROC32_DT<-kwAllPairsDunnTest(PROC32~region, data=metrics_coordi, method='bh')
#Significance at EC-LY, EC-NC
D3264_DT<-kwAllPairsDunnTest(D32_64~region, data=metrics_coordi, method="bh")
#Significance may exist at EC-LY
D12_DT<-kwAllPairsDunnTest(D1_2~region, data=metrics_coordi, method="bh")
#Significance at LY-NC, KT-NC
library(rcompanion)
library(rempsyc)
#Save table
PROC32_DT<-PMCMRTable(PROC32_DT)
D3264_DT<-PMCMRTable(D3264_DT)
D12_DT<-PMCMRTable(D12_DT)
DT_sum<-cbind(D12_DT, D3264_DT, PROC32_DT)
DT_sum[c(3,5)]<-NULL
colnames(DT_sum)[2:4]<-c('p.value of D (1 - 2 cm) cm', 'p.value of D (32 - 64 cm)', 'p.value of PROC - 32 cm')
DT_NT<-nice_table(DT_sum)
DT_NTComparison | p.value of D (1 - 2 cm) cm | p.value of D (32 - 64 cm) | p.value of PROC - 32 cm |
|---|---|---|---|
LY - KT = 0 | 1 | 1 | 0.928 |
LD - KT = 0 | 1 | 0.72 | 1 |
EC - KT = 0 | 1 | 0.563 | 0.928 |
NC - KT = 0 | 0.0537 | 1 | 0.928 |
LD - LY = 0 | 1 | 0.163 | 0.928 |
EC - LY = 0 | 1 | 0.0994 | 0.0411 |
NC - LY = 0 | 0.0127 | 1 | 1 |
EC - LD = 0 | 1 | 1 | 0.75 |
NC - LD = 0 | 0.348 | 0.547 | 0.928 |
NC - EC = 0 | 0.411 | 0.385 | 0.0263 |
flextable::save_as_docx(DT_NT, path = "Figures/PairwiseComparison_indicators_Dunn.docx")
####Plot relationships between regions###
library(ggpubr)
#V4
box_V4_region<-ggplot(data=metrics_coordi,
aes(x=region, y=V4 ))+
geom_boxplot(aes(fill=region))+
labs(y='VRM - 4 cm')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
stat_compare_means(label.y = 0.1, cex=5)
box_V4_region
ggsave('Figures/box_v4_region.jpeg', box_V4_region,
width=15, height=15, unit='cm')
#PROC32
box_PROC32_region<-ggplot(data=metrics_coordi,
aes(x=region, y=PROC32 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label = paste('KW chi-squared= ', round(PROC32_kruskal$statistic, 3)))+
labs(y='PROC - 32 cm')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
box_PROC32_region<-box_PROC32_region+
stat_compare_means(label.y = 0.17, cex=5)+ # Add global p-value
annotate('text', x=1:5, y=0.75, label=c('ns', '*', 'ns', '', '*'),
cex=5)
box_PROC32_region
ggsave('Figures/box_PROC32_region.jpeg', box_PROC32_region,
width=15, height=15, unit='cm')
#PLC4
box_PLC4_region<-ggplot(data=metrics_coordi,
aes(x=region, y=PLC4 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
labs(y='PLC - 4 cm')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
stat_compare_means(label.y = 21, cex=5)
box_PLC4_region
ggsave('Figures/box_PLC4_region.jpeg', box_PLC4_region,
width=15, height=15, unit='cm')
#PLC16
box_PLC16_region<-ggplot(data=metrics_coordi,
aes(x=region, y=PLC16 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(PLC16_kruskal$statistic, 3)))+
labs(y='PLC - 16 cm')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
stat_compare_means(label.y = 5.5, cex=5)
box_PLC16_region
ggsave('Figures/box_PLC16_region.jpeg', box_PLC16_region,
width=15, height=15, unit='cm')
#PLC32
box_PLC32_region<-ggplot(data=metrics_coordi,
aes(x=region, y=PLC32 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(PLC32_kruskal$statistic, 3)))+
labs(y='PLC - 32 cm')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
stat_compare_means(label.y =3.15, cex=4)
box_PLC32_region
ggsave('Figures/box_PLC32_region.jpeg', box_PLC32_region,
width=15, height=15, unit='cm')
#D12
box_D12_region<-ggplot(data=metrics_coordi,
aes(x=region, y=D1_2 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(D12_kruskal$statistic, 3)))+
labs(y='D(1 - 2 cm)')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
box_D12_region<-box_D12_region+
stat_compare_means(label.y = 2.32, cex=5)+ # Add global p-value
annotate('text', x=1:5, y=2.28, label=c('ns', '*', 'ns','ns',''),
cex=5)
box_D12_region
ggsave('Figures/box_D12_region.jpeg', box_D12_region,
width=15, height=15, unit='cm')
#D24
box_D24_region<-ggplot(data=metrics_coordi,
aes(x=region, y=D2_4 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(D24_kruskal$statistic, 3)))+
labs(y='D(2 - 4 cm) cm')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
stat_compare_means(label.y = 2.06, cex=5)
box_D24_region
ggsave('Figures/box_D24_region.jpeg', box_D24_region,
width=15, height=15, unit='cm')
#D816
box_D816_region<-ggplot(data=metrics_coordi,
aes(x=region, y=D8_16 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
labs(y='D(8 - 16 cm)')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
stat_compare_means(label.y = 2.19, cex=5)
box_D816_region
ggsave('Figures/box_D816_region.jpeg', box_D816_region,
width=15, height=15, unit='cm')
#D1632
box_D1632_region<-ggplot(data=metrics_coordi,
aes(x=region, y=D16_32 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
labs(y='D(16 - 32 cm)')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))+
stat_compare_means(label.y = 2.23, cex=5)
box_D1632_region
ggsave('Figures/box_D1632_region.jpeg', box_D1632_region,
width=15, height=15, unit='cm')
#D3264
box_D3264_region<-ggplot(data=metrics_coordi,
aes(x=region, y=D32_64 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
labs(y='D(32 - 64 cm)')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
box_D3264_region<-box_D3264_region+
stat_compare_means(label.y = 2.4, cex=5)+ # Add global p-value
# Pairwise comparison against reference
annotate('text', x=1:5, y=2.36, label=c('ns', 'ns', 'ns', '', 'ns'),
cex=5)
box_D3264_region
ggsave('Figures/box_D3264_region.jpeg', box_D3264_region,
width=15, height=15, unit='cm')
#Sq
box_VRSD_region<-ggplot(data=metrics_coordi,
aes(x=region, y=VRSD4 ))+
geom_boxplot(aes(fill=region))+
ggtitle(label =paste('KW chi-squared= ', round(PLC4_kruskal$statistic, 3)))+
labs(y='Sq')+xlab('Region')+
guides(fill = guide_legend(title = "Region"))+
scale_fill_manual(values=c( "#fc0f00", "#fb7200","#fec700", "#15e0fa","#0c59fe" ))
box_VRSD_region<-box_VRSD_region+
stat_compare_means(label.y = 0.65, cex=5)
box_VRSD_region
ggsave('Figures/box_VRSD_region.jpeg', box_VRSD_region,
width=15, height=15, unit='cm')4.3. PCA based on complexity metrics
PCA was used to visualize the relationships between sites, PERMANOVA was used to examined the interregional differences of complexity.
###################PCA based on metrics#########################
library(vegan)## 載入需要的套件：permute## 載入需要的套件：lattice## This is vegan 2.6-4metr_sel_std<-decostand(metrics_selected, 'standardize')
metrics_pca<-rda(metr_sel_std, scale=T)
PCAsummary<-summary(metrics_pca)
gps<-c('NC', 'NC', 'LY', 'NC', 'LD', 'LD', 'LY', 'KT', 'EC', 'LD', 'LD', 'KT', 'LY',
'EC', 'EC', 'KT', 'EC', 'EC', 'NC', 'KT', 'NC', 'LD', 'LY', 'KT', 'LY')
my_cols<-c("#15e0fa", "#fc0f00","#fec700", "#fb7200" ,"#0c59fe" )
gps2<-as.factor(gps)
levels(gps2)<-my_cols
#Visualization
ordiplot(metrics_pca, type='n',display='site',
xlab=paste('PC1,',
round(PCAsummary[["cont"]][["importance"]][3,1]*100, 2), '%'),
ylab=paste('PC2,',
round(PCAsummary[["cont"]][["importance"]][2,2]*100, 2), '%'),
scaling=2)
pca.scores<-scores(metrics_pca,
choices=1:2,
scaling=2,
display='site')
points(x=pca.scores[,1], y=pca.scores[,2], pch=20, cex=1.5,
col=as.vector(gps2))
ordispider(metrics_pca,
groups=gps,
display='site',
label=F,
lwd=2,
col=c("#15e0fa", "#fc0f00","#fec700", "#fb7200" ,"#0c59fe" ),
scaling=2)
ordihull(metrics_pca,
groups=gps,
display='site',
draw='lines',
lwd = 2,
col=c("#15e0fa", "#fc0f00","#fec700", "#fb7200" ,"#0c59fe" ),
scaling=2)
vct.scores<-scores(metrics_pca,
choices=1:2,
scaling=2,
display='sp')
arrows(x0=0, y0=0,
x1=vct.scores[,1],
y1=vct.scores[,2],
length=0.1,
lty=1,
lwd=3,
col=rgb(0.8, 0.1,0,0.7))
text(x=vct.scores[1,1]+0.3, y= vct.scores[1,2]+0.2,
'VRM - 4 cm', col= 'red', cex=1)
text(x=vct.scores[2,1]-0.3, y= vct.scores[2,2]+0.2,
'PLC - 4 cm', col= 'red', cex=1)
text(x=vct.scores[3,1]+0.5, y= vct.scores[3,2],
'PROC - 32 cm', col= 'red', cex=1)
text(x=vct.scores[4,1]+0.28, y= vct.scores[4,2]+0.26,
'PLC - 32 cm', col= 'red', cex=0.8)
text(x=vct.scores[5,1]-0.8, y= vct.scores[5,2]+0.4,
'PLC - 16 cm', col= 'red', cex=1)
text(x=vct.scores[6,1]+0.3, y= vct.scores[6,2]+0.3,
'D (1 - 2 cm)', col= 'red', cex=1)
text(x=vct.scores[7,1]+0.3, y= vct.scores[7,2]+0.2,
'D (2 - 4 cm)', col= 'red', cex=1)
text(x=vct.scores[8,1]+0.62, y= vct.scores[8,2]+0.1,
'D (8 - 16 cm)', col= 'red', cex=1)
text(x=vct.scores[9,1]+0.75, y= vct.scores[9,2],
'D (16 - 32 cm)', col= 'red', cex=1)
text(x=vct.scores[10,1]+0.55, y= vct.scores[10,2]-0.1,
'D (32 - 64 cm)', col= 'red', cex=1)
text(x=vct.scores[11,1]+0.2, y= vct.scores[11,2]-0.4,
'Sq', col= 'red', cex=1)
legend('topleft', fill=my_cols,
legend= c('East Coast','Kenting', 'Ludao', 'Lanyu', 'North Coast'),
title='Region', cex=1)
#PERMANOVA on complexity
metr_sel_euc<-vegdist(metr_sel_std, 'euclidean')
metr_sel_region<-cbind(metr_sel_std, gps)
metr_sel_PERMANOVA<-adonis2(metr_sel_std~gps, data=metr_sel_region,
permutations = 999, method='euclidean')
metr_sel_PERMANOVA## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = metr_sel_std ~ gps, data = metr_sel_region, permutations = 999, method = "euclidean")
## Df SumOfSqs R2 F Pr(>F)
## gps 4 79.403 0.30077 2.1507 0.007 **
## Residual 20 184.597 0.69923
## Total 24 264.000 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PERMANOVA_nice<-nice_table(metr_sel_PERMANOVA)
flextable::save_as_docx(PERMANOVA_nice, path = "Figures/Complexity_PERMANOVA.docx")
library(RVAideMemoire)## *** Package RVAideMemoire v 0.9-81-2 ***##
## 載入套件：'RVAideMemoire'## 下列物件被遮斷自 'package:raster':
##
## cvlibrary(devtools)## 載入需要的套件：usethis##
## 載入套件：'devtools'## 下列物件被遮斷自 'package:permute':
##
## checklibrary(pairwiseAdonis)## 載入需要的套件：cluster#Pairwise permutation test
permtest<-pairwise.adonis(as.data.frame(metr_sel_std), metr_sel_region$gps,
sim.method='euclidean')
permtest_nice<-nice_table(permtest,note='* p < .05, ** p < .01, *** p < .001', highlight = T)
permtest_nicepairs | Df | SumsOfSqs | F.Model | R2 | p | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
NC vs LY | 1.00 | 18.07 | 1.49 | .16 | .179 | 1.00 | |
NC vs LD | 1.00 | 24.46 | 2.04 | .20 | .050 | 0.50 | |
NC vs KT | 1.00 | 22.74 | 2.01 | .20 | .076 | 0.76 | |
NC vs EC | 1.00 | 40.54 | 3.90 | .33 | .008** | 0.08 | |
LY vs LD | 1.00 | 16.27 | 1.82 | .19 | .120 | 1.00 | |
LY vs KT | 1.00 | 6.09 | 0.74 | .08 | .590 | 1.00 | |
LY vs EC | 1.00 | 35.73 | 4.84 | .38 | .008** | 0.08 | |
LD vs KT | 1.00 | 6.85 | 0.84 | .10 | .485 | 1.00 | |
LD vs EC | 1.00 | 13.04 | 1.81 | .18 | .110 | 1.00 | |
KT vs EC | 1.00 | 14.72 | 2.25 | .22 | .067 | 0.67 | |
Note. * p < .05, ** p < .01, *** p < .001 | |||||||
flextable::save_as_docx(permtest_nice, path = "Figures/Complexity_PairwisePermTest.docx")5. Further simplification on composition data and conducting PCA
Several categories were further merged (i.e., combine Heliopora, Millepora, Tubipora with hard corals). And combining biotic and abiotic (substrates) datasets again. Using PCA to visualize relative relationships between sites, using PERMANOVA to examine the interregional differences.
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
library(corrplot)
library(car)
library(readxl)
library(olsrr)
#Read metric data
StrMetrics<-read.csv('Sites_MetricSummary.csv', header=T)
row.names(StrMetrics)<-StrMetrics[,1]
metrics<-StrMetrics[, 2:24]
#Import relief data
relief_data<-read_xlsx('VRandVRSD.xlsx', sheet=2)
relief_data<-as.data.frame(relief_data)
row.names(relief_data)<-relief_data[,1]
relief_data<-relief_data[-1]
metrics<-merge(metrics, relief_data[2], by='row.names')
row.names(metrics)<-metrics$Row.names
metrics<-metrics[-1]
write.csv(metrics, 'Site_MetricSummary_2.csv')
#########################################################################
########Import composition (Both biotic components and substrates) and further organization###########################
library(dplyr)##
## 載入套件：'dplyr'## 下列物件被遮斷自 'package:raster':
##
## intersect, select, union## 下列物件被遮斷自 'package:car':
##
## recode## 下列物件被遮斷自 'package:stats':
##
## filter, lag## 下列物件被遮斷自 'package:base':
##
## intersect, setdiff, setequal, unioncomposition<-read.csv('analysis_tables/Composition__CombineSpoencAndCru_AllSites.csv',row.names = 1 ) #Biotic components
substrates<-read.csv( 'analysis_tables/Substrates_AllSites.csv') #substrates
substrates<-substrates[,-1]
composition_substrates<-cbind(composition, substrates[,-c(1,6)])
composition_substrates$bss<-1-rowSums(composition_substrates[,3:45]) #get BSS
composition_substrates$hcbus<-composition_substrates$hcbus+composition_substrates$helbus+
composition_substrates$milbus+composition_substrates$hcdig #Get bushy HC
composition_substrates$helbus<-NULL
composition_substrates$milbus<-NULL
composition_substrates$hcdig<-NULL
#get encrusting HC
composition_substrates$hcenc<-composition_substrates$hcenc+
composition_substrates$milenc+composition_substrates$helenc+
composition_substrates$tubenc
composition_substrates$milenc<-NULL
composition_substrates$helenc<-NULL
composition_substrates$tubenc<-NULL
#get massive HC
composition_substrates$hcmas<-composition_substrates$hcmas+composition_substrates$tubmas
composition_substrates$tubmas<-NULL
#get Arborescent HC
composition_substrates$hcarb<-composition_substrates$hcarb+composition_substrates$milarb
composition_substrates$milarb<-NULL
#get macroalgae
composition_substrates$ma<-composition_substrates$com+
composition_substrates$art+composition_substrates$cof
composition_substrates$com<-NULL
composition_substrates$art<-NULL
composition_substrates$cof<-NULL
#get US
composition_substrates$us<-composition_substrates$sandus+
composition_substrates$rubbus+composition_substrates$siltus+
composition_substrates$gravus
composition_substrates$sandus<-NULL
composition_substrates$rubbus<-NULL
composition_substrates$siltus<-NULL
composition_substrates$gravus<-NULL
head(composition_substrates)## site region actinia ceralg hcarb hcbus hccol
## 1 TS KT 1.382217e-03 0.001840512 7.365073e-03 0.07568280 0.001178526
## 2 LDS KT 0.000000e+00 0.000000000 1.601930e-02 0.11848354 0.014356737
## 3 JLS KT 4.594236e-05 0.000148033 5.171541e-02 0.13114650 0.000000000
## 4 HBH KT 0.000000e+00 0.000000000 0.000000e+00 0.05456729 0.007572080
## 5 DiBS KT 0.000000e+00 0.000000000 5.505862e-05 0.07930534 0.002333195
## 6 GW LD 0.000000e+00 0.000000000 0.000000e+00 0.06344037 0.000000000
## hcenc hcfol hcmas hctab ocdig oclob
## 1 0.16715161 0.0050774656 0.06404287 0.040455649 0.0001649510 3.451637e-03
## 2 0.31648276 0.0313315459 0.04485327 0.000000000 0.0124688964 2.069143e-03
## 3 0.03559424 0.0423277998 0.00963379 0.002503439 0.0000000000 1.443082e-03
## 4 0.04454544 0.4555850206 0.01511464 0.000000000 0.0000000000 3.233330e-03
## 5 0.10917749 0.0002641661 0.01333520 0.000000000 0.0003398285 8.729073e-05
## 6 0.18696224 0.2107142169 0.02108280 0.000000000 0.0002013752 0.000000e+00
## ocmas tws fil ocbus ascidian ocfan
## 1 0.0021837048 3.937554e-04 0.0000000000 0.000000000 0.000000000 0.000000e+00
## 2 0.0009843722 0.000000e+00 0.0018783695 0.005894327 0.000000000 0.000000e+00
## 3 0.0000000000 0.000000e+00 0.0035133109 0.002369318 0.001000274 4.835259e-05
## 4 0.0000000000 0.000000e+00 0.0006558537 0.005005998 0.000000000 0.000000e+00
## 5 0.0002757979 9.899104e-05 0.0000000000 0.007079813 0.000000000 0.000000e+00
## 6 0.0000000000 3.013859e-04 0.0000000000 0.000000000 0.000000000 0.000000e+00
## zoaenc spoglo spomas zoamas occlu ocenc
## 1 0.0000000000 0.000000e+00 0.0000000000 0.0000000000 0.0000000 0
## 2 0.0000000000 0.000000e+00 0.0000000000 0.0000000000 0.0000000 0
## 3 0.0002374461 0.000000e+00 0.0000000000 0.0000000000 0.0000000 0
## 4 0.0000000000 0.000000e+00 0.0000000000 0.0000000000 0.0000000 0
## 5 0.0008747468 2.266706e-05 0.0001298605 0.0002267812 0.0000000 0
## 6 0.0000000000 0.000000e+00 0.0000000000 0.0009696243 0.1732429 0
## other sporep spopap coeenc cru_or_spoenc bss ma
## 1 3.117624e-03 0 0 0 0.0004121527 0.4972918 0.0681254510
## 2 8.293761e-05 0 0 0 0.0058026699 0.4038181 0.0254740217
## 3 7.965084e-05 0 0 0 0.0057743596 0.6127960 0.0830581394
## 4 0.000000e+00 0 0 0 0.0009811873 0.1511659 0.0010077042
## 5 1.610879e-03 0 0 0 0.0101201119 0.6875589 0.0693698048
## 6 0.000000e+00 0 0 0 0.0091574471 0.3329709 0.0009566711
## us
## 1 0.06068224
## 2 0.00000000
## 3 0.01656496
## 4 0.26056556
## 5 0.01773404
## 6 0.00000000##Combine metrics with composition
#Import metrics
metrics<-read.csv('Site_MetricSummary_2.csv', row.names = 1, header=T)
row.names(composition_substrates)<-composition_substrates$site
comp_sub_metrics<-merge(composition_substrates, metrics, by='row.names')
row.names(comp_sub_metrics)<-comp_sub_metrics$Row.names
comp_sub_metrics<-comp_sub_metrics[,-1]
#Change colnames
indicators<-c('S - 4 cm', 'TRI - 4 cm', 'VRM - 4 cm',
'S - 32 cm', 'TRI - 32 cm', 'VRM - 32 cm',
'PROC - 4 cm', 'PLC - 4 cm', 'PROC - 32 cm',
'PLC - 32 cm', 'S - 16 cm', 'TRI - 16 cm', 'VRM - 16 cm',
'PROC - 16 cm', 'PLC - 16 cm', 'D(1 - 64 cm)', 'SC',
'D (1 - 2 cm)', 'D (2 - 4 cm)', 'D (4 - 8 cm)',
'D(8 - 16 cm)', 'D(16 - 32 cm)', 'D(32 - 64 cm)',
'Sq')
colnames(comp_sub_metrics)[34:57]<-indicators
#####################PCA for composition#########################################
##############set color###########################
library(vegan)
library(fishualize)
library(devtools)
#devtools::install_github("nschiett/fishualize", force = TRUE)
my_cols<-c("#15e0fa", "#fc0f00","#fec700", "#fb7200" ,"#0c59fe" )
##########################tb-PCA (hellinger transformed)###########################################
#Hellinger transformation on composition
comp_sub_hel<-decostand(composition_substrates[,3:33], method='hellinger')
comp_sub_hel<-decostand(comp_sub_hel, method='standardize')
colnames(comp_sub_hel)<-c('Actinia', 'Ceratodictyon', 'Arborescent HC',
'Bushy HC', 'Columnar HC', 'Encrusting HC',
'Foliose HC', 'Massive HC', 'Tabular HC', 'Digitate OC',
'Lobate OC', 'Massive OC', 'TWS','Filamentous algae', 'Bushy OC',
'Ascidian','Fan OC', 'Encrusting Zoanthids',
'Globular Sponges', 'Massive Sponges','Massive Zoanthids',
'Clustered OC','Encrusting OC', 'Others', 'Repent Sponges',
'Papillate Sponges', 'Corallimorpharia', 'CCA',
'BSS', 'Macroalgae', 'US')
groups<-composition_substrates$region
comp_sub_RDA<-rda(comp_sub_hel)
RDA_summary<-summary(comp_sub_RDA)
#Visualization of PCA
ordiplot(comp_sub_RDA, type='n',
display=c('site', 'sp'),
scaling=3,
xlab=paste('PC1, ',
round(RDA_summary$cont$importance[2,1]*100, 2), '%'),
ylab=paste('PC2, ',round(RDA_summary$cont$importance[2,2]*100, 2), '%'),
)
title(sub='Benthic composition based on sites (hellinger transformed)')
ordihull(comp_sub_RDA,
groups=groups,
draw='lines',
scaling=3,
col=my_cols,
lwd=2,
alpha=0.8)
vct.scores<-scores(comp_sub_RDA,
choices=1:2,
display='sp',
scaling=3)
arrows(x0=0, y0=0,
x1=vct.scores[,1],
y1=vct.scores[,2],
length=0,
lty=1,
lwd=0.7,
col='red')
text(vct.scores[3,1],vct.scores[3,2]-0.05,
row.names(vct.scores)[3],
col='red', cex=0.8) #Arborescent HC
text(vct.scores[4,1]+0.13,vct.scores[4,2],
row.names(vct.scores)[4],
col='red', cex=0.8) #Bushy HC
text(vct.scores[5,1]+0.11,vct.scores[5,2]-0.02,
row.names(vct.scores)[5],
col='red', cex=0.6) #Columnar HC
text(vct.scores[6,1]+0.03,vct.scores[6,2]-0.06,
row.names(vct.scores)[6],
col='red', cex=0.8) #encrusting HC
text(vct.scores[7,1],vct.scores[7,2]-0.02,
row.names(vct.scores)[7],
col='red', cex=0.8) #Foliose HC
text(vct.scores[8,1]-0.12,vct.scores[8,2]-0.05,
row.names(vct.scores)[8],
col='red', cex=0.8) #Massive HC
text(vct.scores[9,1]+0.09,vct.scores[9,2]-0.02,
row.names(vct.scores)[9],
col='red', cex=0.6) #Tabular HC
text(vct.scores[10,1]-0.1,vct.scores[10,2]+0.04,
row.names(vct.scores)[10],
col='red', cex=0.8) #Digitate OC
text(vct.scores[11,1],vct.scores[11,2]-0.02,
row.names(vct.scores)[11],
col='red', cex=0.8) #Lobate OC
text(vct.scores[12,1]+0.1,vct.scores[12,2]-0.05,
row.names(vct.scores)[12],
col='red', cex=0.6) #Massive OC
text(vct.scores[13,1],vct.scores[13,2]+0.02,
row.names(vct.scores)[13],
col='red', cex=0.6) #TWS
text(vct.scores[14,1]+0.03,vct.scores[14,2]+0.04,
row.names(vct.scores)[14],
col='red', cex=0.8)#turf
text(vct.scores[15,1]-0.1,vct.scores[15,2]-0.03,
row.names(vct.scores)[15],
col='red', cex=0.8) #Bushy OC
text(vct.scores[16,1],vct.scores[16,2],
row.names(vct.scores)[16],
col='red', cex=0.6) #can remove
text(vct.scores[18,1]-0.02,vct.scores[18,2]+0.02,
row.names(vct.scores)[18],
col='red', cex=0.8)#zonenc
text(vct.scores[20,1],vct.scores[20,2]+0.02,
row.names(vct.scores)[20],
col='red', cex=0.6) #Massive sponges
text(vct.scores[21,1]+0.02,vct.scores[21,2]+0.02,
row.names(vct.scores)[21],
col='red', cex=0.8) #massive zoanthids
text(vct.scores[22,1]-0.02,vct.scores[22,2]-0.1,
row.names(vct.scores)[22],
col='red', cex=0.55) #clustered oc
text(vct.scores[24,1]-0.05,vct.scores[24,2]+0.05,
row.names(vct.scores)[24],
col='red', cex=0.5)# others
text(vct.scores[26,1]-0.15,vct.scores[26,2]-0.01,
row.names(vct.scores)[26],
col='red', cex=0.6) #papillate sponge
text(vct.scores[28,1]+0.06,vct.scores[28,2]+0.13,
row.names(vct.scores)[28],
col='red', cex=0.8) #CCA
text(vct.scores[29,1],vct.scores[29,2]+0.03,
row.names(vct.scores)[29],
col='red', cex=0.8) #BSS
text(vct.scores[30,1]+0.15,vct.scores[30,2],
row.names(vct.scores)[30],
col='red', cex=0.8) #MA
text(vct.scores[31,1]+0.07,vct.scores[31,2],
row.names(vct.scores)[31],
col='red', cex=0.8) #US
groups2<-groups
groups2<-as.factor(groups2)
levels(groups2)<-my_cols
groups2<-as.vector(groups2)
points(comp_sub_RDA, display='site',
scaling=3,
col=groups2,
pch=20, cex=1.5)
vct.scores1<-scores(comp_sub_RDA,
choices=1:2,
display='site',
scaling=3)
text(vct.scores1[3,1]+0.02,
vct.scores1[3,2]-0.03,
'Jialeshuei',
col="#fc0f00",
cex=0.6
)
text(vct.scores1[20,1]-0.02,
vct.scores1[20,2]+0.04,
'Fenniaolin',
col="#15e0fa",
cex=0.6)
text(vct.scores1[22,1],
vct.scores1[22,2]+0.03,
'Bitou',
col="#0c59fe" ,
cex=0.6)
text(vct.scores1[21,1]+0.08,
vct.scores1[21,2]+0.03,
'Longdong',
col="#0c59fe" ,
cex=0.6)
legend('topleft', fill=my_cols,
legend= c('East Coast','Kenting', 'Ludao', 'Lanyu', 'North Coast'),
title='Region', cex=1.2)
###PERMANOVA on benthic composition##################################
library(rempsyc)
comp_sub_std<-decostand(composition_substrates[,3:33], 'standardize')
comp_sub_hel_std<-decostand(comp_sub_hel, 'standardize')
comp_sub_hel_euc<-vegdist(comp_sub_hel_std, 'euc')
comp_sub_PERMANOVA2<-adonis2(comp_sub_hel_std~region, data=composition_substrates,
permutations = 999, method='euc')
comp_sub_PERMANOVA2## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = comp_sub_hel_std ~ region, data = composition_substrates, permutations = 999, method = "euc")
## Df SumOfSqs R2 F Pr(>F)
## region 4 190.77 0.2564 1.7241 0.001 ***
## Residual 20 553.23 0.7436
## Total 24 744.00 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1comp_sub_PERMANOVA_nice<-nice_table(comp_sub_PERMANOVA2)
flextable::save_as_docx(comp_sub_PERMANOVA_nice, path = "Figures/Composition_PERMANOVA.docx")
#Pairwise Permutation
library(RVAideMemoire)
library(devtools)
#install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
library(pairwiseAdonis)
permtest2<-pairwise.adonis(as.data.frame(comp_sub_hel_std), composition_substrates$region,
sim.method = 'euc')
permtest2<-as.data.frame(permtest2)
permtest2$sig<-NULL
permtest_nice2<-nice_table(permtest2,
note='* p < .05, ** p < .01, *** p < .001', highlight = T)
permtest_nice2pairs | Df | SumsOfSqs | F.Model | R2 | p | p.adjusted |
|---|---|---|---|---|---|---|
KT vs LD | 1.00 | 38.42 | 1.46 | .15 | .085 | 0.85 |
KT vs LY | 1.00 | 35.30 | 1.28 | .14 | .099 | 0.99 |
KT vs EC | 1.00 | 40.99 | 1.18 | .13 | .194 | 1.00 |
KT vs NC | 1.00 | 66.34 | 2.18 | .21 | .009** | 0.09 |
LD vs LY | 1.00 | 32.68 | 1.59 | .17 | .054 | 0.54 |
LD vs EC | 1.00 | 54.13 | 1.95 | .20 | .030* | 0.30 |
LD vs NC | 1.00 | 58.00 | 2.48 | .24 | .006** | 0.06 |
LY vs EC | 1.00 | 54.72 | 1.88 | .19 | .017* | 0.17 |
LY vs NC | 1.00 | 40.37 | 1.63 | .17 | .014* | 0.14 |
EC vs NC | 1.00 | 55.97 | 1.75 | .18 | .009** | 0.09 |
Note. * p < .05, ** p < .01, *** p < .001 | ||||||
flextable::save_as_docx(permtest_nice2, path = "Figures/Composition_PairwisePermTest.docx")6. Initial selection on benthic categories (sparsity + cross-validation)
##################CV by glmnet###########################
comp_sub_metrics<-comp_sub_metrics[,-26] #remove other motile invertebrates
#Set filter of sparsity = 0.5
library(glmnet)## 載入需要的套件：Matrix##
## 載入套件：'Matrix'## 下列物件被遮斷自 'package:tidyr':
##
## expand, pack, unpack## Loaded glmnet 4.1-7filter <- function(x, ...) which(colMeans(x == 0) > 0.5)
sparsity <- function(fraction = 0.7) {
function(x, ...) which(colMeans(x == 0) > fraction)
}
#V4
cv_V4<-cv.glmnet(y=comp_sub_metrics$`VRM - 4 cm`,
x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
exclude=sparsity(0.5),
family=gaussian(),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldV4_coef<-as.data.frame(as.matrix(coef(cv_V4, s='lambda.min')))
#PROC32
cv_PROC32<-cv.glmnet(y=comp_sub_metrics$`PROC - 32 cm`,
x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
exclude=sparsity(0.5),
family=Gamma('log'),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldPROC32_coef<-as.data.frame(as.matrix(coef(cv_PROC32, s='lambda.min')))
#PLC4
cv_PLC4<-cv.glmnet(y=comp_sub_metrics$`PLC - 4 cm`,
x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
exclude=sparsity(0.5),
family=gaussian(),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldPLC4_coef<-as.data.frame(as.matrix(coef(cv_PLC4, s='lambda.min')))
#PLC16
cv_PLC16<-cv.glmnet(y=comp_sub_metrics$`PLC - 16 cm`,
x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
exclude=sparsity(0.5),
family=Gamma('log'),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldPLC16_coef<-as.data.frame(as.matrix(coef(cv_PLC16, s='lambda.min')))
#PLC32
cv_PLC32<-cv.glmnet(y=comp_sub_metrics$`PLC - 32 cm`,
x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
exclude=sparsity(0.5),
family=gaussian(),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldPLC32_coef<-as.data.frame(as.matrix(coef(cv_PLC32, s='lambda.min')))
#D12
cv_D12<-cv.glmnet(y=comp_sub_metrics$`D (1 - 2 cm)`,
x=as.matrix(sqrt(comp_sub_metrics[,3:32])),
exclude=sparsity(0.5),
family=Gamma('log'),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldD12_coef<-as.data.frame(as.matrix(coef(cv_D12, s='lambda.min')))
#D24
cv_D24<-cv.glmnet(y=comp_sub_metrics$`D (2 - 4 cm)`,
x=as.matrix(comp_sub_metrics[,3:32]),
exclude=sparsity(0.5),
family=Gamma('log'),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldD24_coef<-as.data.frame(as.matrix(coef(cv_D24, s='lambda.min')))
#D816
cv_D816<-cv.glmnet(y=comp_sub_metrics$`D(8 - 16 cm)`,
x=as.matrix(comp_sub_metrics[,3:32]),
exclude=sparsity(0.5),
family=Gamma('log'),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldD816_coef<-as.data.frame(as.matrix(coef(cv_D816, s='lambda.min')))
#D1632
cv_D1632<-cv.glmnet(y=comp_sub_metrics$`D(16 - 32 cm)`,
x=as.matrix(comp_sub_metrics[,3:32]),
exclude=sparsity(0.5),
family=gaussian(),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldD1632_coef<-as.data.frame(as.matrix(coef(cv_D1632, s='lambda.min')))
#D3264
cv_D3264<-cv.glmnet(y=comp_sub_metrics$`D(32 - 64 cm)`,
x=as.matrix(comp_sub_metrics[,3:32]),
exclude=sparsity(0.5),
family=gaussian(),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldD3264_coef<-as.data.frame(as.matrix(coef(cv_D3264, s='lambda.min')))
#Sq
cv_VRSD<-cv.glmnet(y=comp_sub_metrics$Sq,
x=as.matrix(comp_sub_metrics[,3:32]),
exclude=sparsity(0.5),
family=gaussian(),
nfolds=25)## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## foldSq_coef<-as.data.frame(as.matrix(coef(cv_VRSD, s='lambda.min')))#Show all plots
par(mfrow=c(2,3))
plot(cv_V4)
title('VRM - 4 cm', line= 3)
plot(cv_PLC4)
title('PLC - 4 cm', line = 3)
plot(cv_D12)
title('D (1 - 2 cm)', line = 3)
plot(cv_D24)
title('D (2 - 4 cm)', line = 3)
plot(cv_PLC16)
title('PLC - 16 cm', line = 3)
plot(cv_D816)
title('D (8 - 16 cm)', line = 3)
par(mfrow=c(2,3))
plot(cv_D1632)
title('D (16 - 32 cm)', line = 3)
plot(cv_PROC32)
title('PROC - 32 cm', line = 3)
plot(cv_PLC32)
title('PLC - 32 cm', line = 3)
plot(cv_D3264)
title('D (32 - 64 cm)', line = 3)
plot(cv_VRSD)
title('Sq', line = 3)
###Show all coefficients
cv_coef<-cbind(V4_coef, PLC4_coef, D12_coef, D24_coef,
PLC16_coef, D816_coef, D1632_coef,
PROC32_coef, PLC32_coef, D3264_coef,
Sq_coef)
colnames(cv_coef)<-c('VRM - 4 cm', 'PLC - 4 cm', 'D (1 - 2 cm)',
'D (2 - 4 cm)', 'PLC - 16 cm', 'D (8 - 16 cm)', 'D (16 - 32 cm)',
'PROC - 32 cm', 'PLC - 32 cm', 'D (32 - 64 cm)', 'Sq' )
row.names(cv_coef)<-c('Intercept', 'Actinia', 'Ceratodictyon', 'Arborescent HC',
'Bushy HC', 'Columnar HC', 'Encrusting HC', 'Foliose HC',
'Massive HC', 'Tabular HC', 'Digitate OC', 'Lobate OC',
'Massive OC', 'TWS', 'Filamentous algae', 'Bushy OC',
'Ascidian', 'Fan OC', 'Encrusting zoanthid', 'Globular sponge',
'Massive sponge', 'Massive zoanthid', 'Clustered OC', 'Encrusting OC',
'Repent sponge', 'Papillate sponge', 'Corallimorpharian', 'CCA',
'BSS', 'Macroalgae', 'US')
cv_coef<-cv_coef[c(1,4,5,7,8,9,11,12,15,16,19,22,28, 29:31 ),]
categories_retained<-row.names(cv_coef)
cv_coef<-cbind(categories_retained,cv_coef)
colnames(cv_coef)[1]<-'Category'
cv_glm_nice<-nice_table(cv_coef)
cv_glm_niceCategory | VRM - 4 cm | PLC - 4 cm | D (1 - 2 cm) | D (2 - 4 cm) | PLC - 16 cm | D (8 - 16 cm) | D (16 - 32 cm) | PROC - 32 cm | PLC - 32 cm | D (32 - 64 cm) | Sq |
|---|---|---|---|---|---|---|---|---|---|---|---|
Intercept | 0.05 | 14.47 | 0.86 | 0.75 | 1.32 | 0.76 | 2.15 | -0.94 | 1.98 | 2.23 | 0.33 |
Arborescent HC | 0.05 | 0.00 | 0.08 | 0.36 | 0.00 | 0.00 | 0.00 | 0.81 | 0.00 | 0.00 | 0.00 |
Bushy HC | 0.05 | 0.00 | 0.05 | 0.10 | 0.00 | 0.02 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
Encrusting HC | 0.00 | 0.00 | -0.11 | -0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
Foliose HC | 0.00 | 0.00 | -0.03 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
Massive HC | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.61 | 0.00 | 0.28 | 0.00 |
Digitate OC | 0.00 | 0.00 | -0.08 | 0.00 | 0.00 | 0.01 | 0.00 | 0.00 | 0.00 | 0.67 | 0.00 |
Lobate OC | 0.05 | 0.00 | 0.04 | 0.00 | 0.00 | 0.00 | 0.00 | 0.43 | 0.00 | 0.00 | 0.00 |
Filamentous algae | 0.03 | 0.00 | -0.04 | 0.00 | 0.00 | 0.00 | 0.00 | -0.57 | 0.00 | -0.06 | 0.00 |
Bushy OC | 0.00 | 0.00 | 0.02 | 0.44 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
Encrusting zoanthid | 0.01 | 0.00 | 0.15 | 0.00 | 0.80 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
Massive zoanthid | 0.00 | 0.00 | -0.08 | 0.00 | 0.16 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
CCA | -0.03 | 0.00 | -0.00 | 0.00 | 0.37 | -0.13 | -0.37 | -1.44 | 0.00 | -0.93 | 0.00 |
BSS | -0.00 | 0.00 | -0.09 | -0.00 | 0.00 | 0.00 | -0.03 | 0.00 | 0.00 | -0.02 | 0.00 |
Macroalgae | 0.00 | 0.00 | -0.04 | 0.00 | 0.00 | 0.00 | 0.00 | 0.07 | 0.00 | 0.00 | 0.00 |
US | 0.00 | 0.00 | -0.04 | 0.00 | 0.00 | 0.00 | 0.00 | 0.67 | 0.00 | 0.07 | 0.22 |
flextable::save_as_docx(cv_glm_nice, path = "Figures/CV_Coefficients.docx")7. Second selection (stepwise VIF) and tb-RDA
Using stepwise VIF selection for the second selection on benthic categories.
Using hellinger-transformed composition to explain standardized complexity. Using permutation test to reveal significance. Using variance partitioning analysis to reveal the importance of each category.
par(mfrow=c(1,1))
library(dplyr)
#Standardize the hellinger-transformed composition
comp_sub_hel_std<-decostand(comp_sub_hel, 'standardize')
#Select the categories retained from initial selection
MetricsRDA_table<-dplyr::select(comp_sub_hel_std,`Arborescent HC`, `Bushy HC`, `Encrusting HC`,`Foliose HC`, `Massive HC`,
`Encrusting Zoanthids`,`Massive Zoanthids`,
`Lobate OC`, `Digitate OC`,`Bushy OC`, `Filamentous algae`, `Macroalgae`,
`CCA`,`BSS`, `US`
)
#Standardize the metrics
metrics_std<-decostand(metrics[,c(3,8,9,10,15,18,19,21,22,23,24)], 'standardize')
rda_table<- merge(MetricsRDA_table, metrics_std, by='row.names')
row.names(rda_table)<-rda_table$Row.names
rda_table<-rda_table[,-1]
rda_mod<-rda(rda_table[,c(16:26)]~., data=rda_table[,1:15])
rda_mod_summary<-summary(rda_mod)
#Second selection on benthic categories (stepwise VIF)
vif.cca(rda_mod)## `Arborescent HC` `Bushy HC` `Encrusting HC`
## 2.759531 2.134464 6.744055
## `Foliose HC` `Massive HC` `Encrusting Zoanthids`
## 8.163998 14.994587 11.260313
## `Massive Zoanthids` `Lobate OC` `Digitate OC`
## 5.858777 5.550216 4.022515
## `Bushy OC` `Filamentous algae` Macroalgae
## 4.606060 3.227758 9.063365
## CCA BSS US
## 4.247677 19.328454 3.120272#remove BSS
#Select the retained categories again
MetricsRDA_table<-dplyr::select(comp_sub_hel_std,`Arborescent HC`, `Bushy HC`, `Encrusting HC`,`Foliose HC`, `Massive HC`,
`Encrusting Zoanthids`,`Massive Zoanthids`,
`Lobate OC`, `Digitate OC`,`Bushy OC`, `Filamentous algae`, `Macroalgae`,
`CCA`, `US`
)
metrics_std<-decostand(metrics[,c(3,8,9,10,15,18,19,21,22,23,24)], 'standardize')
rda_table<- merge(MetricsRDA_table, metrics_std, by='row.names')
row.names(rda_table)<-rda_table$Row.names
#Create a group, showing regions
groups<-c('NC', 'NC', 'LY', 'NC', 'LD', 'LD', 'LY', 'KT', 'EC', 'LD', 'LD',
'KT', 'LY', 'EC', 'EC', 'KT', 'EC', 'EC', 'NC', 'KT', 'NC', 'LD',
'LY', 'KT', 'LY')
groups2<-as.factor(groups)
levels(groups2)<-my_cols
rda_table<-rda_table[,-1]
rda_mod<-rda(rda_table[,c(15:25)]~., data=rda_table[,1:14])
rda_mod_summary<-summary(rda_mod)
vif.cca(rda_mod) #All VIF <10## `Arborescent HC` `Bushy HC` `Encrusting HC`
## 2.599292 1.923347 5.332395
## `Foliose HC` `Massive HC` `Encrusting Zoanthids`
## 2.594163 3.685568 8.328774
## `Massive Zoanthids` `Lobate OC` `Digitate OC`
## 5.541295 4.589677 3.215497
## `Bushy OC` `Filamentous algae` Macroalgae
## 4.006500 3.174003 6.898325
## CCA US
## 2.960635 2.550567##Visualization
#Show R^2
RsquareAdj(rda_mod)## $r.squared
## [1] 0.7542331
##
## $adj.r.squared
## [1] 0.4101595par(mfrow=c(1,1))
ordiplot(rda_mod,type='n',
scaling=3,,
xlab=paste('RDA1, ',
round(rda_mod_summary$concont$importance[2,1]*100, 2), '%'),
ylab=paste('RDA2, ',round(rda_mod_summary$concont$importance[2,2]*100, 2), '%'),
xlim=c(-1,1),
cex=0.8)
title(
sub='Variable selection by removing variables with VIF > 10')
rda_sc<-scores(rda_mod,choice=1:2,
display='sp',
scaling=2)
arrows(x0=0, y0=0,
x1=rda_sc[,1],
y1=rda_sc[,2],
length=0.1,
lty=1,
lwd=1.5,
col='red')
rda_sc1<-scores(rda_mod,choice=1:2,
display='site',
scaling=3)
points(x=rda_sc1[,1], y=rda_sc1[,2], pch=20, cex=1.5,
col=as.vector(groups2))
ordihull(rda_mod,
groups= groups,
draw='lines',
scaling=3,
alpha=0.5,
lwd=1.5,
col=my_cols,
border = NULL)
rda_sc2<-scores(rda_mod,choice=1:2,
scaling=2,
display = 'bp')
row.names(rda_sc2)<-c('Arborescent HC', 'Bushy HC', 'Encrusting HC',
'Foliose HC', 'Massive HC',
'Encrusting Zoanthids','Massive Zoanthids',
'Lobate OC', 'Digitate OC','Bushy OC', 'Filamentous algae',
'Macroalgae','CCA', 'US')
arrows(x0=0, y0=0,
x1=rda_sc2[,1],
y1=rda_sc2[,2],
length=0,
lty=1,
lwd=2,
col='black')
text(x=rda_sc2[1,1]+0.15,
y=rda_sc2[1,2]+0.065,
row.names(rda_sc2)[1],
cex=0.7,
col='black')
text(x=rda_sc2[2,1]+0.265,
y=rda_sc2[2,2]+0.1,
row.names(rda_sc2)[2],
cex=0.8,
col='black')
text(x=rda_sc2[3,1]-0.13,
y=rda_sc2[3,2]-0.12,
row.names(rda_sc2)[3],
cex=0.8,
col='black')
text(x=rda_sc2[4,1]+0.25,
y=rda_sc2[4,2]-0.07,
row.names(rda_sc2)[4],
cex=0.8,
col='black')
text(x=rda_sc2[5,1]+0.18,
y=rda_sc2[5,2]-0.2,
row.names(rda_sc2)[5],
cex=0.8,
col='black')
text(x=rda_sc2[6,1]-0.46,
y=rda_sc2[6,2]-0.017,
row.names(rda_sc2)[6],
cex=0.6,
col='black')
text(x=rda_sc2[6,1]-0.3,
y=rda_sc2[6,2]+0.05,
row.names(rda_sc2)[7],
cex=0.6,
col='black')
text(x=rda_sc2[8,1]+0.37,
y=rda_sc2[8,2]+0.07,
row.names(rda_sc2)[8],
cex=0.8,
col='black')
text(x=rda_sc2[9,1]+0.35,
y=rda_sc2[9,2],
row.names(rda_sc2)[9],
cex=0.8,
col='black')
text(x=rda_sc2[10,1]+0.02,
y=rda_sc2[10,2]+0.1,
row.names(rda_sc2)[10],
cex=0.8,
col='black')
text(x=rda_sc2[11,1]-0.025,
y=rda_sc2[11,2]+0.11,
row.names(rda_sc2)[11],
cex=0.6,
col='black')
text(x=rda_sc2[12,1]+0.11,
y=rda_sc2[12,2]-0.18,
row.names(rda_sc2)[12],
cex=0.55,
col='black')
text(x=rda_sc2[13,1]-0.1,
y=rda_sc2[13,2]+0.05,
row.names(rda_sc2)[13],
cex=0.8,
col='black')
text(x=rda_sc2[14,1]+0.05,
y=rda_sc2[14,2]-0.13,
row.names(rda_sc2)[14],
cex=0.8,
col='black')
text(x=rda_sc[1,1]+0.3,
y=rda_sc[1,2]+0.05,
'VRM - 4 cm',
cex=0.9,
col='red')
text(x=rda_sc[2,1]-0.1,
y=rda_sc[2,2]+0.25,
'PLC - 4 cm',
cex=0.9,
col='red')
text(x=rda_sc[3,1]+0.45,
y=rda_sc[3,2]-0.1,
'PROC - 32 cm',
cex=0.9,
col='red')
text(x=rda_sc[4,1]+0.45,
y=rda_sc[4,2]-0.2,
'PLC - 32 cm',
cex=0.8,
col='red')
text(x=rda_sc[5,1]-0.4,
y=rda_sc[5,2]-0.04,
'PLC - 16 cm',
cex=0.9,
col='red')
text(x=rda_sc[6,1]+0.18,
y=rda_sc[6,2]+0.2,
'D (1 - 2 cm)',
cex=0.9,
col='red')
text(x=rda_sc[7,1]+0.15,
y=rda_sc[7,2]+0.2,
'D (2 - 4 cm)',
cex=0.9,
col='red')
text(x=rda_sc[8,1]+0.37,
y=rda_sc[8,2]+0.05,
'D (8 - 16 cm)',
cex=0.9,
col='red')
text(x=rda_sc[9,1]+0.5,
y=rda_sc[9,2]-0.21,
'D (16 - 32 cm)',
cex=0.9,
col='red')
text(x=rda_sc[10,1]+0.4,
y=rda_sc[10,2]-0.31,
'D (32 - 64 cm)',
cex=0.9,
col='red')
text(x=rda_sc[11,1]+0.15,
y=rda_sc[11,2]-0.1,
'Sq',
cex=0.9,
col='red')
legend('topright', fill=my_cols,
legend= c('East Coast','Kenting', 'Ludao', 'Lanyu', 'North Coast'),
title='Region',
cex=1.2)
#Show summary of the constrained model
rda_mod_summary ##
## Call:
## rda(formula = rda_table[, c(15:25)] ~ `Arborescent HC` + `Bushy HC` + `Encrusting HC` + `Foliose HC` + `Massive HC` + `Encrusting Zoanthids` + `Massive Zoanthids` + `Lobate OC` + `Digitate OC` + `Bushy OC` + `Filamentous algae` + Macroalgae + CCA + US, data = rda_table[, 1:14])
##
## Partitioning of variance:
## Inertia Proportion
## Total 11.000 1.0000
## Constrained 8.297 0.7542
## Unconstrained 2.703 0.2458
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## Eigenvalue 3.3298 1.8034 1.2039 1.00719 0.44074 0.22889 0.105438
## Proportion Explained 0.3027 0.1639 0.1094 0.09156 0.04007 0.02081 0.009585
## Cumulative Proportion 0.3027 0.4666 0.5761 0.66765 0.70772 0.72853 0.738116
## RDA8 RDA9 RDA10 RDA11 PC1 PC2
## Eigenvalue 0.086751 0.049569 0.034350 0.0066187 1.2629 0.6193
## Proportion Explained 0.007886 0.004506 0.003123 0.0006017 0.1148 0.0563
## Cumulative Proportion 0.746002 0.750509 0.753631 0.7542331 0.8690 0.9253
## PC3 PC4 PC5 PC6 PC7 PC8
## Eigenvalue 0.31843 0.21142 0.11673 0.091287 0.035766 0.02728
## Proportion Explained 0.02895 0.01922 0.01061 0.008299 0.003251 0.00248
## Cumulative Proportion 0.95429 0.97351 0.98412 0.992420 0.995672 0.99815
## PC9 PC10
## Eigenvalue 0.011693 0.0086399
## Proportion Explained 0.001063 0.0007854
## Cumulative Proportion 0.999215 1.0000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## Eigenvalue 3.3298 1.8034 1.2039 1.0072 0.44074 0.22889 0.10544
## Proportion Explained 0.4013 0.2174 0.1451 0.1214 0.05312 0.02759 0.01271
## Cumulative Proportion 0.4013 0.6187 0.7638 0.8852 0.93833 0.96592 0.97863
## RDA8 RDA9 RDA10 RDA11
## Eigenvalue 0.08675 0.049569 0.03435 0.0066187
## Proportion Explained 0.01046 0.005975 0.00414 0.0007978
## Cumulative Proportion 0.98909 0.995062 0.99920 1.0000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 4.03089
##
##
## Species scores
##
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## V4 0.9858 0.3700 0.08551 -0.4348 0.048884 -0.13822
## PLC4 -0.4461 0.6915 -0.39389 -0.3200 0.136029 0.27246
## PROC32 0.8833 -0.4579 -0.06117 0.1519 -0.072617 0.12693
## PLC32 0.3694 -0.2888 -0.89377 0.3381 0.049175 0.08735
## PLC16 -0.3837 -0.1067 -0.25535 -0.4762 -0.688945 0.01808
## D1_2 0.7352 0.6647 0.23228 0.5269 -0.187704 -0.07290
## D2_4 0.6346 0.8562 -0.15508 0.1743 -0.224811 0.04522
## D8_16 0.8193 0.1599 -0.10335 -0.5615 0.204370 0.15831
## D16_32 0.7653 -0.4450 -0.29084 -0.3423 0.008021 -0.27649
## D32_64 0.6335 -0.5439 0.07099 0.2073 -0.128538 0.26853
## VRSD4 0.2863 -0.2234 0.75478 -0.2207 -0.060545 0.19929
##
##
## Site scores (weighted sums of species scores)
##
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 82.5 -1.35323 -2.00650 1.930786 -0.43590 -1.82578 0.464597
## BT -0.56245 0.17865 -0.287361 -1.74634 -1.55646 0.855593
## CCG 0.04190 1.73336 2.070078 0.04737 0.74665 0.602983
## CJ 0.07692 -0.73282 0.844272 -0.31376 1.71572 0.890174
## CK 0.63277 -1.32448 0.393205 0.71133 -0.42394 -0.006186
## DaBS -0.90423 -0.37032 -0.489845 0.20567 1.01853 -0.727456
## DC 0.23044 0.89952 -1.198009 -0.32731 -1.50122 -0.071492
## DiBS -0.77649 0.56325 0.574256 1.00468 0.93608 -0.735997
## FNL 1.08432 0.44154 0.005434 -0.12650 0.09923 -0.447854
## GGB 0.22928 -0.29771 -0.238655 0.07425 -0.33688 1.821578
## GW -0.63263 -0.17139 -2.017465 0.89600 0.18060 -0.879072
## HBH 0.45026 0.03335 0.734766 0.89933 -0.50650 0.817287
## HR 0.73786 1.06863 -0.459736 -0.05137 -1.45310 0.278436
## HS 0.92370 -1.54142 -0.617109 -0.07919 0.45972 1.066678
## JH 0.36925 -0.70379 0.622378 -0.17409 1.32887 0.179194
## JLS 1.23257 0.68498 0.693771 0.33762 -0.03229 0.804830
## JMZ 1.59217 -0.70115 -1.046467 -0.84022 0.34913 -0.367680
## JQ 0.06906 -1.01478 0.607659 0.56581 2.01496 0.406128
## LD 0.05324 0.60557 1.113210 -2.57338 0.21174 -1.304019
## LDS -0.08973 0.66957 -1.450940 -0.03823 0.73776 1.473230
## MYS -2.17321 1.11339 -0.749820 0.15501 0.19930 0.688480
## SL 0.95862 0.17935 0.966558 1.71559 -0.81276 -0.918447
## TL -1.04137 0.19287 -0.149560 0.32607 -0.02802 -2.339783
## TS 0.05251 -0.09990 -1.131960 -0.46018 -0.21007 -1.946109
## YY -1.20152 0.60025 -0.719447 0.22775 -1.31129 -0.605095
##
##
## Site constraints (linear combinations of constraining variables)
##
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 82.5 -0.84269 -2.04012 1.060994 0.01211 -1.04809 0.02771
## BT -0.58393 0.23330 -0.247014 -1.79981 -1.66575 0.92760
## CCG 0.03141 1.89197 1.028257 0.40863 0.79207 0.30029
## CJ -0.42795 -0.32892 0.402549 -0.78744 0.32720 -0.02688
## CK -0.24026 -0.66382 -0.737456 0.76690 -0.96858 -0.62998
## DaBS -1.27160 -0.06655 0.399381 0.38701 0.99274 -0.16350
## DC 0.39397 0.81423 -0.588483 0.06418 -0.98809 0.32128
## DiBS -0.36121 0.08902 0.638646 0.57777 0.56759 -0.56014
## FNL 1.00879 0.55986 -0.003252 -0.10787 -0.03265 -0.61270
## GGB 0.29856 -0.30629 -0.056919 -0.29101 -0.39915 2.40790
## GW -0.09683 -0.08946 -1.598913 0.84229 0.09702 -0.59836
## HBH 0.35403 -0.39126 0.645709 1.15982 0.19666 0.99702
## HR 0.59728 0.49384 0.459178 0.10309 -1.49311 0.05654
## HS 0.17828 -1.17547 -0.201842 -0.13487 0.44699 -0.25824
## JH 0.70265 -0.41194 1.345236 -0.34049 0.70598 0.61629
## JLS 1.02630 1.35732 -0.251769 0.12513 -0.21023 0.35993
## JMZ 1.71508 -0.74890 -0.804286 -0.67829 0.79514 -0.39024
## JQ -0.20941 -0.71948 -0.241369 0.47293 1.10643 0.73041
## LD -0.02623 0.43159 0.827087 -2.38006 0.80155 -1.29372
## LDS -0.10623 0.39280 -1.828138 0.06004 0.69795 0.49406
## MYS -1.76685 0.67139 -0.151129 0.17618 0.58388 0.60627
## SL 1.00704 0.05476 1.294795 1.49313 -0.68790 -0.97080
## TL -0.80006 0.28761 -0.337091 -0.06179 -1.10194 -1.41786
## TS 0.69327 -0.87840 -0.829758 -0.36222 0.47264 -0.44032
## YY -1.27339 0.54292 -0.224414 0.29465 0.01167 -0.48256
##
##
## Biplot scores for constraining variables
##
## RDA1 RDA2 RDA3 RDA4 RDA5
## `Arborescent HC` 0.46608 0.32562 -0.1127738 0.184509 -0.398183
## `Bushy HC` 0.46969 0.42827 -0.1144653 0.218628 0.023679
## `Encrusting HC` -0.06931 -0.25548 -0.0480181 0.051118 -0.351759
## `Foliose HC` 0.36195 -0.07125 -0.1436978 0.349544 0.033711
## `Massive HC` 0.32680 -0.28765 0.2218517 -0.067273 -0.471554
## `Encrusting Zoanthids` -0.19488 0.10841 0.1981194 -0.745849 -0.433235
## `Massive Zoanthids` -0.19148 0.08691 -0.0003583 -0.621176 -0.416900
## `Lobate OC` 0.32639 0.10719 -0.0979112 -0.096798 -0.005932
## `Digitate OC` 0.20161 0.04163 -0.2664732 -0.226901 0.003099
## `Bushy OC` 0.27367 0.34742 -0.0588393 0.237920 -0.116191
## `Filamentous algae` -0.20644 0.16162 0.3125838 -0.702297 0.206525
## Macroalgae 0.13039 -0.12046 -0.2565190 -0.177941 0.223476
## CCA -0.63980 0.39215 0.0171568 -0.145532 -0.514588
## US 0.22531 -0.53728 0.2569424 0.004425 0.210593
## RDA6
## `Arborescent HC` 0.179640
## `Bushy HC` -0.018169
## `Encrusting HC` -0.147289
## `Foliose HC` 0.041287
## `Massive HC` -0.076608
## `Encrusting Zoanthids` -0.037089
## `Massive Zoanthids` 0.148523
## `Lobate OC` -0.002256
## `Digitate OC` 0.321167
## `Bushy OC` -0.031960
## `Filamentous algae` -0.177057
## Macroalgae 0.049085
## CCA -0.032177
## US 0.290866#Show the scores of RDA1 and RDA2
rda_scores<-as.data.frame(rda_sc2)
rda_scores<-round(rda_scores, 3)
rda_scores$Groups<-row.names(rda_scores)
rda_scores<-dplyr::select(rda_scores, Groups, RDA1, RDA2)
library(gt)
rda_scores_gt<-gt(rda_scores)
rda_scores_gt| Groups | RDA1 | RDA2 |
|---|---|---|
| Arborescent HC | 0.466 | 0.326 |
| Bushy HC | 0.470 | 0.428 |
| Encrusting HC | -0.069 | -0.255 |
| Foliose HC | 0.362 | -0.071 |
| Massive HC | 0.327 | -0.288 |
| Encrusting Zoanthids | -0.195 | 0.108 |
| Massive Zoanthids | -0.191 | 0.087 |
| Lobate OC | 0.326 | 0.107 |
| Digitate OC | 0.202 | 0.042 |
| Bushy OC | 0.274 | 0.347 |
| Filamentous algae | -0.206 | 0.162 |
| Macroalgae | 0.130 | -0.120 |
| CCA | -0.640 | 0.392 |
| US | 0.225 | -0.537 |
gtsave(rda_scores_gt, 'Figures/RDA_scores.docx')
###############Check significance of RDA model##############
#Check the significance by axis
rda_aov<-anova.cca(rda_mod, by='axis',
step=5000)
rda_aov## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = rda_table[, c(15:25)] ~ `Arborescent HC` + `Bushy HC` + `Encrusting HC` + `Foliose HC` + `Massive HC` + `Encrusting Zoanthids` + `Massive Zoanthids` + `Lobate OC` + `Digitate OC` + `Bushy OC` + `Filamentous algae` + Macroalgae + CCA + US, data = rda_table[, 1:14])
## Df Variance F Pr(>F)
## RDA1 1 3.3298 16.0118 0.001 ***
## RDA2 1 1.8034 8.6720 0.110
## RDA3 1 1.2039 5.7890 0.568
## RDA4 1 1.0072 4.8433 0.628
## RDA5 1 0.4407 2.1194 0.998
## RDA6 1 0.2289 1.1007 1.000
## RDA7 1 0.1054 0.5070 1.000
## RDA8 1 0.0868 0.4172 1.000
## RDA9 1 0.0496 0.2384 1.000
## RDA10 1 0.0344 0.1652 1.000
## RDA11 1 0.0066 0.0318 1.000
## Residual 13 2.7034
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Significance at RDA1
rda_aov<-as.data.frame(rda_aov)
rda_aov$Axis<-row.names(rda_aov)
rda_aov<-dplyr::select(rda_aov, Axis, Df, Variance, `F`, `Pr(>F)`)
rda_aov<-nice_table(rda_aov, highlight = T)
rda_aovAxis | Df | Variance | F | Pr(>F) |
|---|---|---|---|---|
RDA1 | 1.00 | 3.33 | 16.01 | 0.00 |
RDA2 | 1.00 | 1.80 | 8.67 | 0.11 |
RDA3 | 1.00 | 1.20 | 5.79 | 0.57 |
RDA4 | 1.00 | 1.01 | 4.84 | 0.63 |
RDA5 | 1.00 | 0.44 | 2.12 | 1.00 |
RDA6 | 1.00 | 0.23 | 1.10 | 1.00 |
RDA7 | 1.00 | 0.11 | 0.51 | 1.00 |
RDA8 | 1.00 | 0.09 | 0.42 | 1.00 |
RDA9 | 1.00 | 0.05 | 0.24 | 1.00 |
RDA10 | 1.00 | 0.03 | 0.17 | 1.00 |
RDA11 | 1.00 | 0.01 | 0.03 | 1.00 |
Residual | 13.00 | 2.70 |
flextable::save_as_docx(rda_aov, path = "Figures/RDA_ANOVA_ByAxis.docx")
#Check significance of the whole constrained model
rda_aov_all<-anova.cca(rda_mod,
step=5000) #Show overall significance
rda_aov_all## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = rda_table[, c(15:25)] ~ `Arborescent HC` + `Bushy HC` + `Encrusting HC` + `Foliose HC` + `Massive HC` + `Encrusting Zoanthids` + `Massive Zoanthids` + `Lobate OC` + `Digitate OC` + `Bushy OC` + `Filamentous algae` + Macroalgae + CCA + US, data = rda_table[, 1:14])
## Df Variance F Pr(>F)
## Model 14 8.2966 2.1921 0.002 **
## Residual 10 2.7034
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1rda_aov_all<-nice_table(rda_aov_all)
rda_aov_allDf | Variance | F | Pr(>F) |
|---|---|---|---|
14.00 | 8.30 | 2.19 | 0.00 |
10.00 | 2.70 |
flextable::save_as_docx(rda_aov_all, path = "Figures/RDA_ANOVA.docx")
#Significantly tells the variation######################Variance partitioning###################################
varpart_plot<-varpart(rda_table[,c(15:25)],
X=rda_table[,1],
rda_table[,2],
rda_table[,3],
rda_table[,4])
varpart_table<-varpart_plot$part$fract[1:4,]
row.names(varpart_table)<-c('Arborescent HC','Bushy HC', 'Encrusting HC', 'Foliose HC' )
varpart_table<-round(varpart_table, 3)
varpart_plot<-varpart(rda_table[,c(15:25)],
X=rda_table[,5],
rda_table[,6],
rda_table[,7],
rda_table[,8])
varpart_table2<-varpart_plot$part$fract[1:4,]
row.names(varpart_table2)<-c('Massive HC','Encrusting zoanthids', 'Massive zoanthids', 'Lobate OC' )
varpart_table2<-round(varpart_table2, 3)
varpart_plot<-varpart(rda_table[,c(15:25)],
X=rda_table[,9],
rda_table[,10],
rda_table[,11],
rda_table[,12])
varpart_table3<-varpart_plot$part$fract[1:4,]
row.names(varpart_table3)<-c('Digitate OC','Bushy OC', 'Filamentous algae', 'Macroalgae' )
varpart_table3<-round(varpart_table3, 3)
varpart_plot<-varpart(rda_table[,c(15:25)],
rda_table[,13],
rda_table[,14] )
varpart_table4<-varpart_plot$part$fract[1:2,]
row.names(varpart_table4)<-c('CCA', 'US' )
varpart_table4<-round(varpart_table4, 3)
colnames(varpart_table4)[2:3]<-c('R.square', 'Adj.R.square')
varpart_all_group<-rbind(varpart_table, varpart_table2, varpart_table3, varpart_table4)
varpart_all_group$Testable<-NULL
varpart_all_group$Category<-row.names(varpart_all_group)
varpart_all_group<-dplyr::select(varpart_all_group, Category, Df, R.square, Adj.R.square)
library(gt)
varpart_gt<-gt(varpart_all_group)
varpart_gt| Category | Df | R.square | Adj.R.square |
|---|---|---|---|
| Arborescent HC | 1 | 0.096 | 0.056 |
| Bushy HC | 1 | 0.104 | 0.065 |
| Encrusting HC | 1 | 0.023 | -0.020 |
| Foliose HC | 1 | 0.058 | 0.017 |
| Massive HC | 1 | 0.062 | 0.021 |
| Encrusting zoanthids | 1 | 0.077 | 0.036 |
| Massive zoanthids | 1 | 0.056 | 0.015 |
| Lobate OC | 1 | 0.039 | -0.003 |
| Digitate OC | 1 | 0.030 | -0.012 |
| Bushy OC | 1 | 0.052 | 0.011 |
| Filamentous algae | 1 | 0.076 | 0.036 |
| Macroalgae | 1 | 0.022 | -0.020 |
| CCA | 1 | 0.162 | 0.125 |
| US | 1 | 0.077 | 0.037 |
gtsave(varpart_gt, 'Figures/Varpart.docx')8. Displaying scatter plot to examine the relationships between complexity metrics and each category
library(glmmTMB)
library(buildmer)
library(lme4)##
## 載入套件：'lme4'## 下列物件被遮斷自 'package:RVAideMemoire':
##
## dummy## 下列物件被遮斷自 'package:raster':
##
## getDatalibrary(ggplot2)
library(tidyr)
#Select categories with sparsity < 0.5
f<-which(colSums(comp_sub_metrics[,3:32]==0)<0.5*nrow(comp_sub_metrics)) #Select categories with sparsity < 0.5
colnames(comp_sub_metrics)[33:56]<-c('S4', 'T4', 'V4', 'S32', 'T32',
'V32', 'PROC4', 'PLC4', 'PROC32',
'PLC32', 'S16', 'T16', 'V16', 'PROC16',
'PLC16', 'D64', 'SC', 'D1_2', 'D2_4', 'D4_8', 'D8_16',
'D16_32', 'D32_64', 'Sq')
glmm_table1<-comp_sub_metrics[,c(2,f+2, 35)]
glmm_table1$region<-as.factor(glmm_table1$region)
col_names<-c('Region', 'Arborescent HC','Bushy HC', 'Encrusting HC',
'Foliose HC', 'Massive HC', 'Digitate OC',
'Lobate OC', 'Filamentous algae', 'Bushy OC',
'Encrusting zoanthids', 'Massive zoanthids', 'CCA', 'BSS', 'Macroalgae',
'US')
#V4
table<-pivot_longer(data=glmm_table1, cols=-c('V4','region'), names_to = 'MFG', values_to ='cover' )
V4_scater<-ggplot(data=table, aes(y=V4,x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols) +
theme(axis.text=element_text(size=10))
#PROC32
glmm_table2<-comp_sub_metrics[,c(2,f+2, 41)]
glmm_table2[2:15]<-glmm_table2[2:15]
#colnames(glmm_table2)[1:16]<-col_names
table<-pivot_longer(data=glmm_table2, cols=-c('PROC32','region'), names_to = 'MFG', values_to ='cover' )
PROC32_scater<-ggplot(data=table, aes(y=PROC32, x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#PLC4
glmm_table12<-comp_sub_metrics[,c(2,f+2, 40)]
glmm_table12[2:15]<-glmm_table12[2:15]
table<-pivot_longer(data=glmm_table12, cols=-c('PLC4','region'), names_to = 'MFG', values_to ='cover' )
PLC4_scater<-ggplot(data=table, aes(y=PLC4, x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#PLC16
glmm_table4<-comp_sub_metrics[,c(2,f+2, 50-2-1)]
glmm_table4[2:15]<-glmm_table4[2:15]
table<-pivot_longer(data=glmm_table4, cols=-c('PLC16','region'), names_to = 'MFG', values_to ='cover' )
PLC16_scater<-ggplot(data=table, aes(y=PLC16, x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#PLC32
glmm_table3<-comp_sub_metrics[,c(2,f+2, 45-2-1)]
glmm_table3[2:15]<-glmm_table3[2:15]
table<-pivot_longer(data=glmm_table3, cols=-c('PLC32','region'), names_to = 'MFG', values_to ='cover' )
PLC32_scater<-ggplot(data=table, aes(y=PLC32, x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#D12
glmm_table5<-comp_sub_metrics[,c(2,f+2, 53-2-1)]
glmm_table5[2:15]<-glmm_table5[2:15]
table<-pivot_longer(data=glmm_table5, cols=-c('D1_2','region'), names_to = 'MFG', values_to ='cover' )
D12_scater<-ggplot(data=table, aes(y=D1_2, x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#D2_4
glmm_table6<-comp_sub_metrics[,c(2,f+2, 54-2-1)]
glmm_table6[2:15]<-glmm_table6[2:15]
table<-pivot_longer(data=glmm_table6, cols=-c('D2_4','region'), names_to = 'MFG', values_to ='cover' )
D24_scater<-ggplot(data=table, aes(y=`D2_4`, x=cover ))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#D8_16
glmm_table8<-comp_sub_metrics[,c(2,f+2, 56-2-1)]
glmm_table8[2:15]<-glmm_table8[2:15]
table<-pivot_longer(data=glmm_table8, cols=-c('D8_16','region'), names_to = 'MFG', values_to ='cover' )
D816_scater<-ggplot(data=table, aes(y=`D8_16`, x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#D16_32
glmm_table9<-comp_sub_metrics[,c(2,f+2, 57-2-1)]
glmm_table9[2:15]<-glmm_table9[2:15]
table<-pivot_longer(data=glmm_table9, cols=-c('D16_32','region'), names_to = 'MFG', values_to ='cover' )
D1632_scater<-ggplot(data=table, aes(y=`D16_32`, x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#D32_64
glmm_table10<-comp_sub_metrics[,c(2,f+2, 58-2-1)]
glmm_table10[2:15]<-glmm_table10[2:15]
table<-pivot_longer(data=glmm_table10, cols=-c('D32_64','region'), names_to = 'MFG', values_to ='cover' )
D3264_scater<-ggplot(data=table, aes(y=`D32_64`, x=cover))+
geom_point(aes(col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#Sq
glmm_table11<-comp_sub_metrics[,c(2,f+2, 59-2-1)]
glmm_table11[2:15]<-glmm_table11[2:15]
table<-pivot_longer(data=glmm_table11, cols=-c('Sq','region'), names_to = 'MFG', values_to ='cover' )
Sq_scater<-ggplot(data=table, aes(y=Sq, x=cover))+
geom_point(aes( col=region))+
geom_smooth(method='lm')+
facet_wrap(vars(MFG), scales='free')+
scale_color_manual(values=my_cols)+
theme(axis.text=element_text(size=10))
#Show all plots
V4_scater## `geom_smooth()` using formula = 'y ~ x'
PROC32_scater## `geom_smooth()` using formula = 'y ~ x'
PLC4_scater## `geom_smooth()` using formula = 'y ~ x'
PLC16_scater## `geom_smooth()` using formula = 'y ~ x'
PLC32_scater## `geom_smooth()` using formula = 'y ~ x'
D12_scater## `geom_smooth()` using formula = 'y ~ x'
D24_scater## `geom_smooth()` using formula = 'y ~ x'
D816_scater## `geom_smooth()` using formula = 'y ~ x'
D1632_scater## `geom_smooth()` using formula = 'y ~ x'
D3264_scater## `geom_smooth()` using formula = 'y ~ x'
Sq_scater## `geom_smooth()` using formula = 'y ~ x'
9.1. Confirmation of relationships by using GLMM
The coefficients were visualized by dwplots.
########################################GLMM########################################
#################################vif prior to modelling######################################
library(dplyr)
library(usdm)
glmm_vif_table<-dplyr::select(glmm_table11,hcarb, hcbus,hcenc,hcfol,
hcmas, zoaenc,zoamas, ocdig,oclob, ocbus,ma, fil, cru_or_spoenc, us)
vifstep(glmm_vif_table) #oclob and zoaenc are removed
###################GLMM modelling###################################################
#####################################################################
colnames(glmm_table1)[17]<-'V4'
V4_glmm<-buildglmmTMB(V4~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table1, family=gaussian(),
buildmerControl = list(include = ~ (1|region), direction='forward'))
V4_mod_sum<-summary(V4_glmm)
colnames(glmm_table12)[c(1,17)]<-c('region','PLC4')
PLC4_glmm<-buildglmmTMB(PLC4~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table12, family=gaussian(),
buildmerControl = list(include = ~ (1|region), direction='forward'))
PLC4_mod_sum<-summary(PLC4_glmm)
#plot(PLC4_glmm)
colnames(glmm_table4)[17]<-'PLC16'
PLC16_glmm<-buildglmmTMB(PLC16~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table4, family=Gamma('log'),
buildmerControl = list(include = ~ (1|region), direction='forward'))
PLC16_mod_sum<-summary(PLC16_glmm)
#plot(PLC32_glmm)
colnames(glmm_table3)[17]<-'PLC32'
PLC32_glmm<-buildglmmTMB(PLC32~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table3, family=gaussian(),
buildmerControl = list(include = ~ (1|region), direction='forward'))
PLC32_mod_sum<-summary(PLC32_glmm)
#plot(PLC32_glmm)
colnames(glmm_table5)[17]<-'D1_2'
D12_glmm<-buildglmmTMB(D1_2~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table5, family=Gamma('log'),
buildmerControl = list(include = ~ (1|region), direction='forward'))
D12_mod_sum<-summary(D12_glmm)
#plot(D12_glmm)
colnames(glmm_table6)[17]<-'D2_4'
D24_glmm<-buildglmmTMB(D2_4~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table6, family=Gamma('log'),
buildmerControl = list(include = ~ (1|region), direction='forward'))
D24_mod_sum<-summary(D24_glmm)
colnames(glmm_table8)[17]<-'D8_16'
D816_glmm<-buildglmmTMB(D8_16~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table8, family=Gamma('log'),
buildmerControl = list(include = ~ (1|region), direction='forward'))
D816_mod_sum<-summary(D816_glmm)
colnames(glmm_table9)[17]<-'D16_32'
D1632_glmm<-buildglmmTMB(D16_32~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table9, family=gaussian(),
buildmerControl = list(include = ~ (1|region), direction='forward'))
D1632_mod_sum<-summary(D1632_glmm)
#plot(D1632_glmm)
colnames(glmm_table10)[17]<-'D32_64'
D3264_glmm<-buildglmmTMB(D32_64~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table10, family=gaussian(),
buildmerControl = list(include = ~ (1|region), direction='forward'))
D3264_mod_sum<-summary(D3264_glmm)
colnames(glmm_table11)[17]<-'Sq'
VRSD_glmm<-buildglmmTMB(Sq~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table11, family=gaussian(),
buildmerControl = list(include = ~ (1|region), direction='forward'))
VRSD_mod_sum<-summary(VRSD_glmm) # No benthic variables can significantly explain Sq
colnames(glmm_table2)[17]<-'PROC32'
PROC32_glmm<-buildglmmTMB(PROC32~hcarb+hcbus+hcenc+hcfol+hcmas+zoamas+ocbus+ocdig+
fil+cru_or_spoenc+ma+us+(1|region),
data=glmm_table2, family=Gamma('log'),
buildmerControl = list(include = ~ (1|region)))
PROC32_mod_sum<-summary(PROC32_glmm)###########Show coefficients#########################
library(broom.mixed)
library(dotwhisker)
V4_mod<-glmmTMB(V4_mod_sum$call$formula,data=glmm_table1,family=gaussian())
V4_mod_sum<-summary(V4_mod)
V4_coef<-dwplot(V4_mod, effects='fixed',whisker_args = list(color = "blue"),
dot_args = list(color = "black") )+
geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('hcbus'='Bushy HC',
'fil'='Filamentous algae',
'ocbus' = 'Bushy OC',
'cru_or_spoenc' = 'CCA',
'hcarb' = 'Arborescent HC',
'zoamas' = 'Massive zoanthids'))+
ggtitle('VRM - 4 cm')
V4_coef
PROC32_mod<-glmmTMB(PROC32_mod_sum$call$formula,data=glmm_table2,family=Gamma('log'))
PROC32_mod_sum<-summary(PROC32_mod)
PROC32_coef<-dwplot(PROC32_mod, effects='fixed',whisker_args = list(color = "red"),
dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA',
'hcmas' = 'Massive HC',
'fil' = 'Filamentous algae'))+
ggtitle('PROC - 32 cm')
PROC32_coef
PLC4_mod<-glmmTMB(PLC4_mod_sum$call$formula,data=glmm_table12,family=gaussian())
PLC4_mod_sum<-summary(PLC4_mod)
PLC4_coef<-dwplot(PLC4_mod, effects='fixed',whisker_args = list(color = "blue"),
dot_args = list(color = "black") )+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels= c('hcmas'='Massive HC','cru_or_spoenc'='CCA'))+
ggtitle('PLC - 4 cm')
PLC4_coef
PLC16_mod<-glmmTMB(PLC16_mod_sum$call$formula,data=glmm_table4,family=Gamma('log'))
PLC16_mod_sum<-summary(PLC16_mod)
PLC16_coef<-dwplot(PLC16_mod, effects='fixed',whisker_args = list(color = "orange"),
dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA',
'ma' = 'Macroalgae'))+
ggtitle('PLC - 16 cm')
PLC16_coef
PLC32_mod<-glmmTMB(PLC32_mod_sum$call$formula,data=glmm_table3,family=gaussian())
PLC32_mod_sum<-summary(PLC32_mod)
PLC32_coef<-dwplot(PLC32_mod, effects='fixed',whisker_args = list(color = "red"),
dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('fil' = 'Filamentous algae'))+
ggtitle('PLC - 32 cm')
PLC32_coef
D12_mod<-glmmTMB(D12_mod_sum$call$formula,data=glmm_table5,family=Gamma('log'))
D12_mod_sum<-summary(D12_mod)
D12_coef<-dwplot(D12_mod, effects='fixed',whisker_args = list(color = "blue"),
dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('hcarb' = 'Arborescent HC',
'hcbus' = 'Bushy HC',
'ocbus' = 'Bushy OC',
'hcfol' = 'Foliose HC',
'hcenc' = 'Encrusting HC',
'us' = 'US',
'zoamas' = 'Massive zoanthids'))+
ggtitle('D (1 - 2 cm)')
D12_coef
D24_mod<-glmmTMB(D24_mod_sum$call$formula,data=glmm_table6,family=Gamma('log'))
D24_mod_sum<-summary(D24_mod)
D24_coef<-dwplot(D24_mod, effects='fixed',whisker_args = list(color = "blue"),
dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('hcarb' = 'Arborescent HC',
'hcbus' = 'Bushy HC',
'ocbus' = 'Bushy OC',
'hcfol' = 'Foliose HC',
'fil' = 'Filamentous algae',
'cru_or_spoenc' = 'CCA'))+
ggtitle('D (2 - 4 cm)')
D24_coef
D816_mod<-glmmTMB(D816_mod_sum$call$formula,data=glmm_table8,family=Gamma('log'))
D816_mod_sum<-summary(D816_mod)
D816_coef<-dwplot(D816_mod, effects='fixed',whisker_args = list(color = "orange"),
dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA',
'zoamas' = 'Massive zoanthids',
'hcfol' = 'Foliose HC',
'fil' = 'Filamentous algae',
'ocdig' = 'Digitate OC'))+
ggtitle('D (8 - 16 cm)')
D816_coef
D1632_mod<-glmmTMB(D1632_mod_sum$call$formula,data=glmm_table9,family=gaussian())
D1632_mod_sum<-summary(D1632_mod)
D1632_coef<-dwplot(D1632_mod, effects='fixed',whisker_args = list(color = "orange"),
dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA'))+
ggtitle('D (16 - 32 cm)')
D1632_coef
D3264_mod<-glmmTMB(D3264_mod_sum$call$formula,data=glmm_table10,family=gaussian())
D3264_mod_sum<-summary(D3264_mod)
D3264_coef<-dwplot(D3264_mod, effects='fixed',whisker_args = list(color = "red"),
dot_args = list(color = "black"))+geom_vline(xintercept=0,lty=2)+
scale_y_discrete(labels=c('cru_or_spoenc' = 'CCA',
'hcmas' = 'Massive HC'))+
ggtitle('D (32 - 64 cm)')
D3264_coef
VRSD_mod<-glmmTMB(VRSD_mod_sum$call$formula,data=glmm_table11,family=gaussian())9.2. Showing the mean value and 95% CI of each coefficient in GLMMs
####################Show CI and coefficients of every model#################################
#V4
V4_CI<-confint(V4_mod)[-nrow(confint(V4_mod)),]
stats_table <- as.data.frame(V4_mod_sum$coefficients$cond)
stats_table_V4<-cbind(
row.names(stats_table),
round(stats_table, 3), V4_CI[,-3]
)
colnames(stats_table_V4) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_V4$Term<-c('Intercept','Bushy HC', 'Filamentous algae',
'Bushy OC', 'CCA', 'Arborescent HC', 'Massive Zoanthids' )
V4_nt<-nice_table(stats_table_V4, highlight = TRUE,
note="* p < .05, ** p < .01, *** p < .001")
V4_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 0.050 | 0.00 | 12.36 | < .001*** | [0.04, 0.06] |
Bushy HC | 0.193 | 0.04 | 5.33 | < .001*** | [0.12, 0.26] |
Filamentous algae | 0.103 | 0.02 | 4.98 | < .001*** | [0.06, 0.14] |
Bushy OC | 0.673 | 0.19 | 3.60 | < .001*** | [0.31, 1.04] |
CCA | -0.219 | 0.07 | -3.20 | .001** | [-0.35, -0.09] |
Arborescent HC | 0.345 | 0.15 | 2.36 | .018* | [0.06, 0.63] |
Massive Zoanthids | 0.075 | 0.03 | 2.22 | .026* | [0.01, 0.14] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#PROC32
PROC32_CI<-confint(PROC32_mod)[-nrow(confint(PROC32_mod)),]
stats_table <- as.data.frame(PROC32_mod_sum$coefficients$cond)
stats_table_PROC32<-cbind(
row.names(stats_table),
round(stats_table,3), PROC32_CI[,-3]
)
colnames(stats_table_PROC32) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_PROC32$Term<-c('Intercept','CCA', 'Massive HC',
'Filamentous algae' )
PROC32_nt<-nice_table(stats_table_PROC32, highlight = TRUE,
note="* p < .05, ** p < .01, *** p < .001")
PROC32_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | -0.755 | 0.10 | -7.40 | < .001*** | [-0.96, -0.56] |
CCA | -7.450 | 1.96 | -3.80 | < .001*** | [-11.30, -3.60] |
Massive HC | 2.008 | 0.81 | 2.48 | .013* | [0.42, 3.59] |
Filamentous algae | -1.411 | 0.50 | -2.80 | .005** | [-2.40, -0.42] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#PLC4
PLC4_CI<-confint(PLC4_mod)[-nrow(confint(PLC4_mod)),]
stats_table <- as.data.frame(PLC4_mod_sum$coefficients$cond)
stats_table_PLC4<-cbind(
row.names(stats_table),
round(stats_table,3), PLC4_CI[,-3]
)
colnames(stats_table_PLC4) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_PLC4$Term<-c('Intercept','Massive HC',
'CCA' )
PLC4_nt<-nice_table(stats_table_PLC4, highlight = TRUE,
note="* p < .05, ** p < .01, *** p < .001")
PLC4_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 14.469 | 0.75 | 19.39 | < .001*** | [13.01, 15.93] |
Massive HC | -19.127 | 6.40 | -2.99 | .003** | [-31.68, -6.57] |
CCA | 46.577 | 15.81 | 2.95 | .003** | [15.60, 77.56] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#PLC16
PLC16_CI<-confint(PLC16_mod)[-nrow(confint(PLC16_mod)),]
stats_table <- as.data.frame(PLC16_mod_sum$coefficients$cond)
stats_table_PLC16<-cbind(
row.names(stats_table),
round(stats_table,3), PLC16_CI[,-3]
)
colnames(stats_table_PLC16) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_PLC16$Term<-c('Intercept','CCA',
'Macroalgae' )
PLC16_nt<-nice_table(stats_table_PLC16, highlight = TRUE,
note="* p < .05, ** p < .01, *** p < .001")
PLC16_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 1.305 | 0.05 | 26.62 | < .001*** | [1.21, 1.40] |
CCA | 3.324 | 1.11 | 3.00 | .003** | [1.15, 5.50] |
Macroalgae | 0.277 | 0.48 | 0.58 | .565 | [-0.66, 1.22] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#PLC32
PLC32_CI<-confint(PLC32_mod)[-nrow(confint(PLC32_mod)),]
stats_table <- as.data.frame(PLC32_mod_sum$coefficients$cond)
stats_table_PLC32<-cbind(
row.names(stats_table),
round(stats_table,3), PLC32_CI[,-3]
)
colnames(stats_table_PLC32) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_PLC32$Term<-c('Intercept','Filamentous algae' )
PLC32_nt<-nice_table(stats_table_PLC32, highlight=T, note="* p < .05, ** p < .01, *** p < .001")
PLC32_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 2.056 | 0.09 | 22.68 | < .001*** | [1.88, 2.23] |
Filamentous algae | -2.633 | 0.94 | -2.79 | .005** | [-4.48, -0.78] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#D12
D12_CI<-confint(D12_mod)[-nrow(confint(D12_mod)),]
stats_table <- as.data.frame(D12_mod_sum$coefficients$cond)
stats_table_D12<-cbind(
row.names(stats_table),
round(stats_table,3), D12_CI[,-3]
)
colnames(stats_table_D12) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_D12$Term<-c('Intercept','Arborescent HC','Bushy HC', 'Bushy OC',
'Foliose HC', 'Encrusting HC', 'US', 'Massive zoanthids')
D12_nt<-nice_table(stats_table_D12, highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D12_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 0.753 | 0.01 | 108.26 | < .001*** | [0.74, 0.77] |
Arborescent HC | 0.566 | 0.19 | 2.98 | .003** | [0.19, 0.94] |
Bushy HC | 0.134 | 0.04 | 3.03 | .002** | [0.05, 0.22] |
Bushy OC | 0.484 | 0.23 | 2.07 | .039* | [0.02, 0.94] |
Foliose HC | 0.058 | 0.03 | 1.89 | .059 | [-0.00, 0.12] |
Encrusting HC | -0.039 | 0.02 | -1.75 | .080 | [-0.08, 0.00] |
US | -0.053 | 0.05 | -1.09 | .274 | [-0.15, 0.04] |
Massive zoanthids | -0.021 | 0.04 | -0.50 | .616 | [-0.10, 0.06] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#D24
D24_CI<-confint(D24_mod)[-nrow(confint(D24_mod)),]
stats_table <- as.data.frame(D24_mod_sum$coefficients$cond)
stats_table_D24<-cbind(
row.names(stats_table),
round(stats_table,3), D24_CI[,-3]
)
colnames(stats_table_D24) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_D24$Term<-c('Intercept','Arborescent HC','Bushy HC', 'Bushy OC',
'Foliose HC', 'Filamentous algae', 'CCA')
D24_nt<-nice_table(stats_table_D24, highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D24_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 0.738 | 0.00 | 223.45 | < .001*** | [0.73, 0.74] |
Arborescent HC | 0.520 | 0.11 | 4.63 | < .001*** | [0.30, 0.74] |
Bushy HC | 0.152 | 0.03 | 5.58 | < .001*** | [0.10, 0.20] |
Bushy OC | 0.721 | 0.14 | 5.09 | < .001*** | [0.44, 1.00] |
Foliose HC | 0.051 | 0.01 | 3.27 | .001** | [0.02, 0.08] |
Filamentous algae | 0.041 | 0.02 | 2.59 | .010* | [0.01, 0.07] |
CCA | 0.110 | 0.05 | 2.12 | .034* | [0.01, 0.21] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#D816
D816_CI<-confint(D816_mod)[-nrow(confint(D816_mod)),]
stats_table <- as.data.frame(D816_mod_sum$coefficients$cond)
stats_table_D816<-cbind(
row.names(stats_table),
round(stats_table,3), D816_CI[,-3]
)
colnames(stats_table_D816) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_D816$Term<-c('Intercept','CCA','Foliose HC', 'Massive zoanthids',
'Filamentous algae', 'Digitate OC')
D816_nt<-nice_table(stats_table_D816, highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D816_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 0.764 | 0.00 | 145.70 | < .001*** | [0.75, 0.77] |
CCA | -0.409 | 0.12 | -3.46 | .001** | [-0.64, -0.18] |
Foliose HC | -0.052 | 0.03 | -1.53 | .125 | [-0.12, 0.01] |
Massive zoanthids | 0.091 | 0.06 | 1.63 | .104 | [-0.02, 0.20] |
Filamentous algae | 0.042 | 0.03 | 1.26 | .207 | [-0.02, 0.11] |
Digitate OC | 0.470 | 0.33 | 1.43 | .154 | [-0.18, 1.12] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#D1632
D1632_CI<-confint(D1632_mod)[-nrow(confint(D1632_mod)),]
stats_table <- as.data.frame(D1632_mod_sum$coefficients$cond)
stats_table_D1632<-cbind(
row.names(stats_table),
round(stats_table,3), D1632_CI[,-3]
)
colnames(stats_table_D1632) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_D1632$Term<-c('Intercept','CCA')
D1632_nt<-nice_table(stats_table_D1632, highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D1632_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 2.147 | 0.01 | 167.63 | < .001*** | [2.12, 2.17] |
CCA | -0.832 | 0.33 | -2.51 | .012* | [-1.48, -0.18] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#D3264
D3264_CI<-confint(D3264_mod)[-nrow(confint(D3264_mod)),]
stats_table <- as.data.frame(D3264_mod_sum$coefficients$cond)
stats_table_D3264<-cbind(
row.names(stats_table),
round(stats_table,3), D3264_CI[,-3]
)
colnames(stats_table_D3264) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_D3264$Term<-c('Intercept', 'CCA', 'Massive HC')
D3264_nt<-nice_table(stats_table_D3264, highlight=T, note="* p < .05, ** p < .01, *** p < .001")
D3264_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 2.224 | 0.02 | 122.72 | < .001*** | [2.19, 2.26] |
CCA | -1.343 | 0.38 | -3.50 | < .001*** | [-2.10, -0.59] |
Massive HC | 0.439 | 0.16 | 2.82 | .005** | [0.13, 0.74] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#Sq
VRSD_CI<-confint(VRSD_mod)[-nrow(confint(VRSD_mod)),]
stats_table <- as.data.frame(VRSD_mod_sum$coefficients$cond)
VRSD_CI<-as.data.frame(VRSD_CI)
stats_table_VRSD<-cbind(
row.names(stats_table),
round(stats_table,3), t(VRSD_CI[-3,])
)
colnames(stats_table_VRSD) <- c(
"Term", "B", "SE", "t", "p",
"CI_lower", "CI_upper")
stats_table_VRSD$Term<-c('Intercept')
VRSD_nt<-nice_table(stats_table_VRSD, highlight=T, note="* p < .05, ** p < .01, *** p < .001")
VRSD_ntTerm | β | SE | t | p | 95% CI |
|---|---|---|---|---|---|
Intercept | 0.343 | 0.03 | 11.23 | < .001*** | [0.28, 0.40] |
Note. * p < .05, ** p < .01, *** p < .001 | |||||
#Save tables
flextable::save_as_docx(V4_nt, path = "Figures/V4_GLMM.docx")
flextable::save_as_docx(PROC32_nt, path = "Figures/PROC32_GLMM.docx")
flextable::save_as_docx(PLC4_nt, path = "Figures/PLC4_GLMM.docx")
flextable::save_as_docx(PLC16_nt, path = "Figures/PLC16_GLMM.docx")
flextable::save_as_docx(PLC32_nt, path = "Figures/PLC32_GLMM.docx")
flextable::save_as_docx(D12_nt, path = "Figures/D12_GLMM.docx")
flextable::save_as_docx(D24_nt, path = "Figures/D24_GLMM.docx")
flextable::save_as_docx(D816_nt, path = "Figures/D816_GLMM.docx")
flextable::save_as_docx(D1632_nt, path = "Figures/D1632_GLMM.docx")
flextable::save_as_docx(D3264_nt, path = "Figures/D3264_GLMM.docx")
flextable::save_as_docx(VRSD_nt, path = "Figures/Sq_GLMM.docx")9.3. Examining the GLMMs (pseudo-Rsquare and model diagnistic)
############################Check pseudo R2##########################################
library(MuMIn)
pseudo_r2_table<-as.data.frame(matrix( ncol=2, nrow=11))
row.names(pseudo_r2_table)<-c('VRM - 4 cm',
'PROC - 32 cm',
'PLC - 4 cm',
'PLC - 16 cm',
'PLC - 32 cm',
'D(1 - 2 cm)',
'D(2 - 4 cm)',
'D(8 - 16 cm)',
'D(16 - 32 cm)',
'D(32 - 64 cm)',
'Sq')
pseudo_r2_table[1,]<-round(r.squaredGLMM(V4_mod), 4) #model with best conditional R2
pseudo_r2_table[2,]<-round(r.squaredGLMM(PROC32_mod)[1,], 4)
pseudo_r2_table[3,]<-round(r.squaredGLMM(PLC4_mod),4)
pseudo_r2_table[4,]<-round(r.squaredGLMM(PLC16_mod)[1,], 4)
pseudo_r2_table[5,]<-round(r.squaredGLMM(PLC32_mod), 4)#The least explained model
pseudo_r2_table[6,]<-round(r.squaredGLMM(D12_mod)[1,],4)
pseudo_r2_table[7,]<-round(r.squaredGLMM(D24_mod)[1,],4)# The second best-explained model
pseudo_r2_table[8,]<-round(r.squaredGLMM(D816_mod)[1,],4)
pseudo_r2_table[9,]<-round(r.squaredGLMM(D1632_mod),4)
pseudo_r2_table[10,]<-round(r.squaredGLMM(D3264_mod),4)
pseudo_r2_table[11,]<-round(r.squaredGLMM(VRSD_mod),4)
colnames(pseudo_r2_table)<-c('Marginal R squared', 'Constrained R squared')
#Save tables
library(gt)
library(rempsyc)
pseudo_r2_table$Indicator<-row.names(pseudo_r2_table)
pseudo_r2_table<-dplyr::select(pseudo_r2_table, Indicator, `Marginal R squared`, `Constrained R squared`)
pseudo_r2_nice<-gt(pseudo_r2_table)
pseudo_r2_nice| Indicator | Marginal R squared | Constrained R squared |
|---|---|---|
| VRM - 4 cm | 0.7648 | 0.7648 |
| PROC - 32 cm | 0.6110 | 0.6626 |
| PLC - 4 cm | 0.4387 | 0.4387 |
| PLC - 16 cm | 0.2758 | 0.2758 |
| PLC - 32 cm | 0.2447 | 0.2447 |
| D(1 - 2 cm) | 0.6515 | 0.6515 |
| D(2 - 4 cm) | 0.8058 | 0.8058 |
| D(8 - 16 cm) | 0.4387 | 0.4387 |
| D(16 - 32 cm) | 0.2281 | 0.2573 |
| D(32 - 64 cm) | 0.4723 | 0.4723 |
| Sq | 0.0000 | 0.0000 |
gtsave(pseudo_r2_nice, 'Figures/PseudoR2.docx')
######################Model diagniostics################################
library(DHARMa)
library(SuppDists)
V4_resi<-simulateResiduals(V4_mod)
plot(V4_resi, sub='VRM - 4 cm model diagnosis')
PROC32_resi<-simulateResiduals(PROC32_mod)
plot(PROC32_resi, sub='PROC - 32 cm model diagnosis')
PLC4_resi<-simulateResiduals(PLC4_mod)
plot(PLC4_resi, sub='PLC - 4 cm model diagnosis')
PLC16_resi<-simulateResiduals(PLC16_mod)
plot(PLC16_resi, sub='PLC - 16 cm model diagnosis')
PLC32_resi<-simulateResiduals(PLC32_mod)
plot(PLC32_resi, sub='PLC - 32 cm model diagnosis')
D12_resi<-simulateResiduals(D12_mod)
plot(D12_resi, sub='D (1 - 2 cm) model diagnosis')
D24_resi<-simulateResiduals(D24_mod)
plot(D24_resi, sub='D (2 - 4 cm) model diagnosis')
D816_resi<-simulateResiduals(D816_mod)
plot(D816_resi, sub='D (8 - 16 cm) model diagnosis')
D1632_resi<-simulateResiduals(D1632_mod)
plot(D1632_resi, sub='D (16 - 32 cm) model diagnosis')
D3264_resi<-simulateResiduals(D3264_mod)
plot(D3264_resi, sub='D (32 - 64 cm) model diagnosis')
VRSD_resi<-simulateResiduals(VRSD_mod)
plot(VRSD_resi, sub='Sq model diagnosis')
10. 3D model quality
Showing the reprojection error, ground resolution, and control error of each model.
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
########################################################################
library(readxl)
mod_qual<-read_xlsx('ModelQuality.xlsx')## New names:
## • `` -> `...1`mod_qual<-as.data.frame(mod_qual)
row.names(mod_qual)<-mod_qual[,1]
mod_qual<-mod_qual[,-1]
mod_qual$Reproj_Error_mm<-mod_qual$Res_mm*mod_qual$Reproj_Error
#Showing resolutions, reprojection error (both in pixel and mm), and planar area of each quadrat.
head(mod_qual)## images Res_mm Reproj_Error Cont_Error PLA Reproj_Error_mm
## 82.5 1545 0.858 0.652 0.009873100 23.9740 0.559416
## BT 1621 0.736 0.554 0.000184904 23.0720 0.407744
## CK 1443 0.998 0.531 0.000607361 24.5164 0.529938
## CJ 1473 0.811 0.591 0.000196119 22.6973 0.479301
## CCG 1657 0.679 0.624 0.000299013 23.6267 0.423696
## DaBS 1527 0.812 0.548 0.000218082 23.7451 0.444976#Showing the ranges of parameters
mod_qual_sum<-data.frame()
for(i in 1:ncol(mod_qual)){
mod_qual_sum[1,i]<-min(mod_qual[,i])+(max(mod_qual[,i])-min(mod_qual[,i]))/2
mod_qual_sum[2,i]<-(max(mod_qual[,i])-min(mod_qual[,i]))/2
mod_qual_sum[3,i]<-mean(mod_qual[,i])
}
colnames(mod_qual_sum)<-colnames(mod_qual)
row.names(mod_qual_sum)<-c('Intermediate', 'HalfRange', 'Mean')
#Showing the ranges
round(mod_qual_sum, 3)## images Res_mm Reproj_Error Cont_Error PLA Reproj_Error_mm
## Intermediate 1394.5 0.875 0.604 0.005 23.792 0.584
## HalfRange 440.5 0.284 0.136 0.005 2.316 0.262
## Mean 1504.6 0.774 0.606 0.001 24.154 0.47311. Displaying the maps
rm(list=ls())
setwd("~/Data_for_Lab_Morris_Wu/myExperiment")
library(maptools)
library(rgdal)
library (sp)
library(raster)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggspatial)
library(rgbif)
#library(mapr)
library(marmap)
library(leaflet)
library(ggsn)
library(RgoogleMaps)
library(ggmap)
library(mapproj)
#install.packages("leaflet")
library(leaflet)
library(devtools)
#Import geographical data of quadrats
library(readxl)
coordi<-read_xlsx('Photogrammetric_Plot_Coordinates.xlsx')
coordi<-as.data.frame(coordi)
coordi<-coordi[,-c(4,5,6)]
row.names(coordi)<-coordi[,1]
#Separate the data by regions
KT<-coordi[c(6:10),]
LY<-coordi[c(11:15),]
LD<-coordi[c(1:5),]
EC<-coordi[c(16:20), ]
NC<-coordi[c(21:25),]
###############################################
library(ggmap)
library(ggsn)
#Installing the key of google map
register_google(key = "AIzaSyCccvogSM3f_4VjI7k7OmvWzA2lAKDUwQc")
#Get map of Taiwan
map <- get_map(location = 'Taiwan', zoom = 7, map='satellite')
#Plot the sampling sites on map
taiwan<-ggmap(map)+
scale_x_continuous(limits = c(119.2, 122.5), expand = c(0, 0)) +
scale_y_continuous(limits = c(21.5, 25.6), expand = c(0, 0))+
scalebar(x.min = 119.5, x.max=120, y.min=22, y.max=22.5,
dist = 30, model = 'WGS84',dist_unit = 'km',
transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
box.fill = c('white', 'yellow'),
st.color='white',
st.bottom = F) +
labs(x='Longitudes', y='Latitudes')+
geom_rect(aes(xmin = 120.6, xmax = 121, ymin = 21.8, ymax = 22.2),
color = '#fc0f00',
fill = NA,lwd=0.8)+
geom_rect(aes(xmin = 121.49, xmax = 121.64, ymin = 21.95, ymax = 22.14),
color = '#F28522',
fill = NA,lwd=0.8)+
geom_rect(aes(xmin = 121.44, xmax = 121.53, ymin = 22.62, ymax = 22.69),
color = "#fec700",
fill = NA,lwd=0.8)+
geom_rect(aes(xmin = 121.22, xmax = 121.85, ymin = 22.86, ymax = 24.49),
color = "#15e0fa",
fill = NA,lwd=0.8)+
geom_rect(aes(xmin = 121.756, xmax = 121.965, ymin = 25, ymax = 25.164),
color = "#0c59fe",
fill = NA, lwd=0.8)+
annotate('text', x=120.5, y=21.7, label='Kenting',size=5,
color='white')+
annotate('text', x=121.5, y=21.9, label='Lanyu',size=5,
color='white')+
annotate('text', x=121.5, y=22.5, label='Ludao',size=5,
color='white')+
annotate('text', x=122, y=22.75, label='East Coast',size=5,
color='white')+
annotate('text', x=122.1, y=25.28, label='North Coast',size=5,
color='white')
taiwan
ggsave('Figures/Taiwan_map.jpg', width=1500, height=2250, units='px')
#Plot map of Kenting
map1 <- get_map(location = c(lon = 120.769022, lat = 21.953408),
zoom = 12, map='satellite')
KT_map<-ggmap(map1)+
scale_x_continuous(limits = c(120.68, 120.88), expand = c(0, 0))+
scale_y_continuous(limits = c(21.88, 22.05), expand = c(0, 0))+
scalebar(x.min = 120.7, x.max=120.75, y.min=21.89, y.max=21.92,
dist = 2, model = 'WGS84',dist_unit = 'km',
transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
box.fill = c('white', 'yellow'),
st.color='white',
st.bottom = F) +
labs(x='Longitudes', y='Latitudes')+
geom_point(data=KT, aes(x=Long, y=Lat), color='#fc0f00',
size=4)+
annotate('text', x=120.71, y=21.937, label='Dinbaisha',size=3,
color='white')+
annotate('text', x=120.755, y=21.923, label='Leidashih',size=3,
color='white')+
annotate('text', x=120.76, y=21.937, label='Houbihu',size=3,
color='white')+
annotate('text', x=120.78, y=21.96, label='Tiaoshih',size=3,
color='white')+
annotate('text', x=120.85, y=21.985, label='Jialeshuei',size=3,
color='white')
KT_map
ggsave('Figures/Kenting_map.jpg', width=1500, height=1500, units='px')
#Plot map of Lanyu
map2 <- get_map(location = c(lon = 121.551093, lat = 22.047084),
zoom = 12, map='satellite')
LY_map<-ggmap(map2)+
scale_x_continuous(limits = c(121.47, 121.625), expand = c(0, 0))+
scale_y_continuous(limits = c(21.99, 22.12), expand = c(0, 0))+
scalebar(x.min = 121.45, x.max=121.5, y.min=22, y.max=22.025,
dist = 1, model = 'WGS84',dist_unit = 'km',
transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
box.fill = c('white', 'yellow'),
st.color='white',
st.bottom = F) +
labs(x='Longitudes', y='Latitudes')+
geom_point(data=LY, aes(x=Long, y=Lat), color='#F28522',
size=4)+
annotate('text', x=121.55, y=22.01, label='Green Grassland',size=3,
color='white')+
annotate('text', x=121.50, y=22.05, label='Yayo',size=3,
color='white')+
annotate('text', x=121.55, y=22.088, label='Hen Rock',size=3,
color='white')+
annotate('text', x=121.575, y=22.094, label='Two Lions',size=3,
color='white')+
annotate('text', x=121.585, y=22.058, label='Donching',size=3,
color='white')
LY_map
ggsave('Figures/Lanyu_map.jpg', width=1500, height=1500, units='px')
#Map of Ludao
map3 <- get_map(location = c(lon = 121.490517 , lat = 22.659331),
zoom = 13, map='satellite')
LD_map<-ggmap(map3)+
scale_x_continuous(limits = c(121.455, 121.52), expand = c(0, 0))+
scale_y_continuous(limits = c(22.625, 22.688), expand = c(0, 0))+
scalebar(x.min = 121.46, x.max=121.48, y.min=22.63, y.max=22.64,
dist = 1, model = 'WGS84',dist_unit = 'km',
transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
box.fill = c('white', 'yellow'),
st.color='white',
st.bottom = F) +
labs(x='Longitudes', y='Latitudes')+
geom_point(data=LD, aes(x=Long, y=Lat), color="#fec700",
size=4)+
annotate('text', x=121.488, y=22.637, label='Dabaisha',size=3,
color='white')+
annotate('text', x=121.48, y=22.642, label='Guiwan',size=3,
color='white')+
annotate('text', x=121.469, y=22.653, label='Shihlang',size=3,
color='white')+
annotate('text', x=121.48, y=22.68, label='Chaikou',size=3,
color='white')+
annotate('text', x=121.493, y=22.683, label='Gongguanbi',size=3,
color='white')
LD_map
ggsave('Figures/Ludao_map.jpg', width=1500, height=1500, units='px')
##Map of East Coast
map4 <- get_map(location = c(lon = 121.608442, lat = 23.934540),
zoom = 8, map='satellite')
EC_map<-ggmap(map4)+
scale_x_continuous(limits = c(121.177, 121.932), expand = c(0, 0))+
scale_y_continuous(limits = c(22.8, 24.58), expand = c(0, 0))+
scalebar(x.min = 121.6, x.max=121.8, y.min=22.9, y.max=23.2,
dist = 10, model = 'WGS84',dist_unit = 'km',
transform = T,st.size = 2.5,st.dist=0.3,box.color = 'white',
box.fill = c('white', 'yellow'),
st.color='white',
st.bottom = F) +
labs(x='Longitudes', y='Latitudes')+
geom_point(data=EC, aes(x=Long, y=Lat), color="#15e0fa",
size=4)+
annotate('text', x=121.35, y=22.88, label='Jiamuziwan',size=3,
color='white')+
annotate('text', x=121.5, y=23.1, label='Jihui',size=3,
color='white')+
annotate('text', x=121.65, y=23.62, label='Hsinshe',size=3,
color='white')+
annotate('text', x=121.62, y=23.7, label='Jiqi',size=3,
color='white')+
annotate('text', x=121.7, y=24.5, label='Fenniaolin',size=3,
color='white')
EC_map
ggsave('Figures/EastCoast_map.jpg', width=1000, height=3000, units='px')
#Map of North Coast
map5 <- get_map(location = c(lon = 121.883732 , lat = 25.112880 ),
zoom = 12, map='satellite')
NC_map<-ggmap(map5)+
scale_x_continuous(limits = c(121.78, 121.95), expand = c(0, 0))+
scale_y_continuous(limits = c(25.05, 25.21), expand = c(0, 0))+
scalebar(x.min = 121.9, x.max=121.945, y.min=25.2, y.max=25.21,
dist = 2, model = 'WGS84',dist_unit = 'km',
transform = T,st.size = 2.5,st.dist=0.35,box.color = 'white',
box.fill = c('white', 'yellow'),
st.color='white',
st.bottom = F) +
labs(x='Longitudes', y='Latitudes')+
geom_point(data=NC, aes(x=Long, y=Lat), color="#0c59fe",
size=4)+
annotate('text', x=121.81, y=25.15, label='Chaojing',size=3,
color='white')+
annotate('text', x=121.895, y=25.128, label='82.5K',size=3,
color='white')+
annotate('text', x=121.91, y=25.14, label='Bitou',size=3,
color='white')+
annotate('text', x=121.92, y=25.11, label='Longdong',size=3,
color='white')+
annotate('text', x=121.9, y=25.075, label='Meiyenshan',size=3,
color='white')
NC_map
ggsave('Figures/NorthCoast_map.jpg', width=1500, height=1500, units='px')
#################################################################